Next Article in Journal
The Migration and Invasion of Oral Squamous Carcinoma Cells: Matrix, Growth Factor and Signalling Involvement
Previous Article in Journal
Circumventing Drug Treatment? Intrinsic Lethal Effects of Polyethyleneimine (PEI)-Functionalized Nanoparticles on Glioblastoma Cells Cultured in Stem Cell Conditions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Data Driven Mathematical Model of FOLFIRI Treatment for Colon Cancer

1
Department of Mathematics and Statistics, University of Massachusetts Amherst, Amherst, MA 01003, USA
2
Department of Mathematics, Tufts University, Medford, MA 02155, USA
*
Author to whom correspondence should be addressed.
Cancers 2021, 13(11), 2632; https://doi.org/10.3390/cancers13112632
Submission received: 31 March 2021 / Revised: 21 May 2021 / Accepted: 21 May 2021 / Published: 27 May 2021
(This article belongs to the Section Cancer Therapy)

Abstract

:

Simple Summary

Since the micro-environment of colonic tumors, including their immune structure would affect the response to treatments, we study the response of five groups of tumors clustered based on their immune patterns to a common colon cancer treatment. We develop a data driven mathematical model to investigate the behavior of key players in colonic tumors in each of these clusters in response to the FOLFIRI treatment. Although the model shows clear differences in the behavior of tumors in different clusters, it cannot suggest a unique optimal treatment strategy for each cluster. The results show that there is not much difference in the dynamics of tumors in response to 5-FU alone versus 5-FU plus Leucovorin. However, adding Irinotecan changes the dynamics of T-reg and dendritic cells leading to a remarkably slower tumor recurrence, especially for tumors in a cluster, which has the highest level of T-reg/T-helper ratio compared to the other clusters.

Abstract

Many colon cancer patients show resistance to their treatments. Therefore, it is important to consider unique characteristic of each tumor to find the best treatment options for each patient. In this study, we develop a data driven mathematical model for interaction between the tumor microenvironment and FOLFIRI drug agents in colon cancer. Patients are divided into five distinct clusters based on their estimated immune cell fractions obtained from their primary tumors’ gene expression data. We then analyze the effects of drugs on cancer cells and immune cells in each group, and we observe different responses to the FOLFIRI drugs between patients in different immune groups. For instance, patients in cluster 3 with the highest T-reg/T-helper ratio respond better to the FOLFIRI treatment, while patients in cluster 2 with the lowest T-reg/T-helper ratio resist the treatment. Moreover, we use ROC curve to validate the model using the tumor status of the patients at their follow up, and the model predicts well for the earlier follow up days.

1. Introduction

Colorectal cancer (CRC), the third most common cancer diagnosed in both men and women in the United States excluding skin cancers, is estimated to cause about 52,980 deaths during 2021 [1]. Surgical resection, radiation therapy and systemic therapies that use medications such as chemotherapy, targeted therapy, and immunotherapy are the treatment options depending on several factors such as the type and stage of the disease, the molecular analysis of the tumor, possible side effects and overall health of the patients [2,3]. Early stage tumors can be curable with surgical resection while many patients with advance stage and metastatic CRC receive chemotherapy as a combination of treatment [4]. Although overall mortality rate of CRC patients has been decreasing for several decades, reduction rate slowed over the past decade (2008–2017) [5]. Survival rate remains poor for patients with metastatic CRC despite advances in the primary treatment of chemotherapy [6]. Predicting variability in response to treatments to increase survival rate and arrive at precision medicine, we need to understand disease progression and determine the major drivers for each patient.
One of the major players in response to cancer treatments is immune system. Immune cells in the tumor microenvironment contact with tumor cells directly or through chemokine and cytokine signaling and they play essential roles in improvement and blocking of therapeutic efficacy and the behavior of the tumor [7]. Targeting tumor cells by radiotherapy and chemotherapy causes the release of damage-associated molecular pattern (DAMP) molecules such as high mobility group box 1 (HMGB1) as a result of necrotic cell death [8,9,10,11], and it has been found that HMGB1 triggers immune responses [12,13]. Activated CD8 + T-cells release a high level of cytokines such as IFN- γ and FasL that boost production of necrotic cells in colon cancer [14]. Dendritic cells are activated by HMGB1 that can be released from macrophages [10] and activated dendritic cells, and they cause activation of T-cells [15]. Moreover, tumor-associated macrophages (TAMs) are known as key regulators of therapeutic response in the tumor microenvironment. CD4 + T-cells release IFN- γ that activates M 1 macrophages [16,17], and CD4 + T-cells are activated by TNF- α , which is released by M 1 macrophages [18]. In contrast, M2 macrophages are activated by exposure to certain cytokines such as IL-4, IL-10 and IL-13 and elevate tumorigenesis [19,20].
Many studies have reported a relationship between clinical outcome and immune cells in colon cancer. For example, a high proportion of CD8+ T-cells, effector memory T-cells and CD4+ T-cells is correlated with longer survival in colorectal cancer [21,22,23]. Furthermore, it has been observed that radiotherapy mediates tumor regression because of the release of IFN- γ by CD8 + T-cells [24,25]. In addition, patients with high expression levels of the Th17 markers show a poor prognosis, while a high expression of the Th1 markers is associated with prolonged disease-free survival in colorectal cancer [26]. Moreover, it has been shown that TAMs mediate resistance to some chemotherapeutic agents such as 5-fluorouracil, doxorubicin [27].
Most chemotherapy treatments for colon cancer include Fluorouracil (5-FU), which is a fluoropyrimidine antimetabolite drug used for different cancers types such as colorectal, breast, head and neck cancers [28]. However, response rate of 5-FU-based chemotherapy as a first-line treatment for advanced colorectal cancer remains only 10–15% [29,30]. To overcome this therapeutic resistance, combinations of chemotherapy drugs such as FOLFIRI (Folinic acid, Fluorouracil and Irinotecan) are used for targeting tumor cells and the tumor microenvirinment simultaneously, and they have improved the response rates up to 40–50% for advanced colorectal cancer [29,31]. Accumulating evidence indicates that the relative abundance of various immune cells and their interaction network with treatment approaches are essential in the colonic tumors’ initiation and progression. Therefore, this study focuses on the interaction between tumor-infiltrating immune cells and FOLFIRI agents by dividing patients into similar cohorts based on their tumor-infiltrating immune variations to model the response to the cancer treatment.
Mathematical models are used in many studies to understand tumor growth dynamics, improve therapeutic responses, find the best treatment combination and overcome drug resistance in different cancer treatments [32,33,34,35,36,37,38,39,40]. The effect of radiotherapy and chemotherapy on tumor growth has been modelled using partial differential equations (PDEs), ordinary differential equations (ODEs) and linear quadratic models in breast and brain tumors [41,42,43]. Immune cell interactions with tumor cells are used as an alternative approach for the mathematical modeling of the cancer treatments in a few studies that pharmacokinetic ODEs are defined to predict the optimal dosing regimen and the combined effect of chemotherapy and immunotherapy [44,45,46,47]. Many of these models such as [45] can not be verified because of lack of time course data for the growth of treated and/or untreated tumors. For this reason, some models such as [46] have simulated outcomes for groups of virtual patients on treatment protocols for which clinical trial data are available, using a range of biologically reasonable patient-specific parameter values.
We have recently developed a data driven mathematical model of colon cancer with a focus on the key players and the interaction network between immune cells and cancer cells in order to discover differences in tumor growths of patients with different immune profiles [48]. In this study, we found that there are five distinct groups of primary colon tumors based on their immune profiles, which have been estimated from their gene expression profiles. To analyze the model’s predictions for tumors in each cluster, patient-specific parameters have been estimated using the data of each cluster [48]. In this study, we extend our previous model including the interaction between Fluorouracil, Leucovorin, Irinotecan, and various cell types in tumor to investigate the effect of these drugs on tumors in each cluster.

2. Materials & Methods

2.1. Mathematical Model

There is a complex web of numerous interactions in the colon cancer microenvironment. To be able to analyse and study the role of key interactions in tumors’ progression, the network of these interactions has been reduced to a clear and compact model in [48], highlighting the key components. This network model contains 14 variables that can be majorly grouped into T-cells, dendritic cells, macrophages, cytokines, cancer cells and necrotic cells. In this paper, we add the interactions of three FOLFIRI agents to this network to study the effect of these drugs on colonic tumors (Figure 1 and Table 1).
The interaction network among some of the key players in colon cancer as modeled in [48] is summarized below. In the model described in [48], individually modelled cytokines include HMGB1 (H), IFN- γ ( I γ ), and TGF- β ( G β ). The group of carcinogenic cytokines including IL-6, IL-17, IL-21 and IL-22 is modeled as one variable μ 1 , while the combined sum of immunosuppressive molecules, including IL-10 and CCL20 is modeled as μ 2 . μ 1 is collectively modelled as secreted by macrophages [49,50,51,52], helper T-cells [51,52,53,54] and a sub-population of dendritic cells [55,56], and μ 2 as produced by macrophages [57,58,59], dendritic cells [55,60] and T-reg cells [53,61,62,63]. HMGB1 is modeled as passively released from necrotic cells [64], or actively secreted from activated T-cells, and macrophages [65,66]. IFN- γ is secreted by a sub-population of macrophages [57,67,68,69,70], helper T-cells [16,17] and cytotoxic cells [14]. TGF- β is produced by macrophages [57,58] and T-reg cells [53,61,62,71].
As it has been described in [48], cells that have been modeled are T-cells, macrophages, dendritic cells, necrotic cells, and cancer cells. Naive T-cells ( T N ) are included in the system to make the system more stable by modeling the activation rates of sub-types of T-cells proportional to the density of naive T-cells. Helper T-cells ( T h ) are modeled as they are activated by dendritic cells or certain cytokines including IL-12, IL-6, and IL-23 [61]. Cytotoxic cells ( T C ) are activated by IL-12, IL-4, IL-5 and IL-13 [15,72,73], and regulatory T-cells ( T r ) are also shown to be activated by the cytokines IL-2 [61,74,75], CCL20 [59], and TGF- β [61,71]. Additionally, T-reg cells inhibit both T h and T C cells by various means, including the production of immunosuppresive cytokines [61]. The major effects on dendritic cells are by HMGB1, which activates [10] them but also reduce their maturation rate as shown by some sources [63,76]; and cancer cells, which indirectly may promote their apoptosis [76,77,78,79,80]. Dendritic cells are modelled as two types of naive ( D N ) and activated (D). Macrophages are either activated by IFN- γ or by interleukins (ILs) [57,81,82]. Proliferation in cancer is taken to be proportional to [ C ] ( 1 [ C ] / C 0 ) , where C 0 [83,84] is the total capacity, with additional proliferation by cytokines [61,85], and IL-6 may cause an additional loss of apoptosis in cancer cells [77,79,86,87]. While, cancer development is suppressed by TGF- β , IL-12, IFN- γ and cytotoxic T-cells [61,88,89,90,91]. Here, necrotic cells are considered to be all those cells that go through the process of necrotic cell death and are modelled with a rate of production given by a fraction of dying cancer cells. This is because they are produced either naturally by the tumor, or due to the indirect effect of cytotoxic T-cells [14].
The ODE system obtained as a result of these interactions (presented in [48]) is:
d T N d t = A T N λ T h D D + λ T h M M + λ T h μ 1 μ 1 T N λ T C T h T h + λ T C D D T N λ T r T h T h + λ T r μ 2 μ 2 + λ T r G β G β T N δ T N T N
d T h d t = λ T h D D + λ T h M M + λ T h μ 1 μ 1 T N δ T h μ 2 μ 2 + δ T h T r T r + δ T h T h
d T C d t = λ T C T h T h + λ T C D D T N δ T C μ 2 μ 2 + δ T C T r T r + δ T C T C
d T r d t = λ T r T h T h + λ T r μ 2 μ 2 + λ T r G β G β T N δ T r μ 1 μ 1 + δ T r T r
d D N d t = A D N λ D H H + λ D C C D N δ D H H + δ D D N
d D d t = λ D H H + λ D C C D N δ D H H + δ D C C + δ D D
d M d t = λ M μ 2 μ 2 + λ M I γ I γ + λ M T h T h M 0 M δ M M
d C d t = λ C + λ C μ 1 μ 1 C 1 C C 0 δ C G β G β + δ C I γ I γ + δ C T C T C + δ C C
d N d t = α N C δ C G β G β + δ C I γ I γ + δ C T C T C + δ C C δ N N
d H d t = λ H N N + λ H M M + λ H T h T h + λ H T C T C + λ H T r T r δ H H
d μ 1 d t = λ μ 1 T h T h + λ μ 1 M M + λ μ 1 D D δ μ 1 μ 1
d μ 2 d t = λ μ 2 M M + λ μ 2 D D + λ μ 2 T r T r δ μ 2 μ 2
d I γ d t = λ I γ T h T h + λ I γ T C T C + λ I γ M M δ I γ I γ
d G β d t = λ G β M M + λ G β T r T r δ G β G β
As it has been mentioned in [48], this system includes 14 variables and 59 parameters, where λ parameters correspond to proliferation, activation and production rates, while δ parameters denote the degradation and cell death rates. The parameters A T N and A D N respectively are the production rates of naive T-cells and dendritic cells (D), and M 0 and C 0 are the total density of macrophages (naive and activated together) and cancer cells’ maximum capacity, respectively.
To this given system, we add the individual interactions between the three drugs of the FOLFIRI regimen - Leucovorin (Folinic Acid), Fluorouracil (5-FU), and Irinotecan; and the variables of the above system. The metabolism and action pathways of these three drugs are complex, involving a number of different molecules and enzymes. Therefore, even though we initially attempted to have a comprehensive system with all of their pharmacodynamical reactions, due to the lack of available parameter values in research, we decided to simplify the model. Our condensed system focuses on the overall change in the drug concentrations and their effect on cells in the tumor microenvironment. The final condensed interaction network between FOLFIRI and the colon cancer environment is given in Figure 1.
Since we are adding the drugs to the system already developed in [48], we employ the same formula as given by the mass action law for defining the terms in our ordinary differential equations. Namely, for any biochemical process A + B C , the differential equation for C is given by d C d t = λ A B , where λ is the production rate of C [48,92,93]. Similarly, an inhibition process D E is given by d E d t = δ D E , where δ is the inhibition rate of E. Additionally, alphas ( α ) are included as constant parameters in the drug equations to necessarily differentiate between the inhibition rates on the targeted cells and the decay effect on the drug concentration itself. In the following, we explain the derived equations for the dynamics of each drug and the cells effected by them, with concentrations given in time, per unit day (changes to the equations from [48] have been highlighted in bold).

2.1.1. Cancer & Necrotic Cells

We model the effects of FOLFIRI drugs (5-FU, Leucovorin, and Irinotecan) on cancer by modifing the Equation (8) in the following way.
d [ C ] d t = ( λ C + λ C μ 1 [ μ 1 ] ) [ C ] ( 1 [ C ] C 0 ) ( δ C G β [ G β ] + δ C I γ [ I γ ] + δ C T C [ T C ] + δ C + δ C 5 fu [ 5 FU ] + δ C 5 fu I γ [ 5 FU ] [ I γ ] δ 5 fu M [ 5 FU ] [ M ] + δ C LV 5 fu [ 5 FU ] [ LV ] + δ C Ir [ Ir ] ) [ C ] ,
where the decay rates represent the following: δ C 5 f u and δ C I r are the direct cytotoxic effects of 5-FU and Irinotecan respectively; δ C 5 f u I γ is the rate of cancer cell death due to the increased activation of 5-FU by IFN- γ ; and lastly δ C L V 5 f u is the combined inhibitory effect of 5-FU and Leucovorin.
Note, 5-FU causes damage to cancer cells by inhibiting essential processes in DNA and RNA synthesis [28,29]. This is done by two main pathways, either by the missincorporation of fluoronucleotides in both RNA and DNA, or by inhibiting the enzyme thymidylate synthase (TYMS), which is a crucial component of DNA replication and repair [28,29,94]. These pathways are aided by IFN- γ , which up-regulates the activities of 5-FU anabolic enzymes, increasing its cytotoxic effect [29]. Hence, Equation (15) includes the terms δ C 5 fu [ 5 FU ] [ C ] , δ C 5 fu I γ [ 5 FU ] [ I γ ] [ C ] . Meanwhile, macrophages have been shown to decrease the inhibitory effect of 5-FU on cancer cells [95,96,97]; modeled by the term δ 5 fu M [ 5 FU ] [ M ] [ C ] in the above Equation (15).
Leucovorin is often administered simultaneously with 5-FU in the treatment of colon cancer, because it has been shown to increase the efficacy of 5-FU and the overall survival rate in patients [94,98,99]. This effect is due to the fact that it helps to further stabilize the complex formed between TYMS and the 5-FU derivative, and thus increasing the retention of 5-FU toxicity [29,100]. Therefore, the effect of Leucovorin on cancer is modeled by the term δ C LV 5 fu [ 5 FU ] [ LV ] [ C ] in Equation (15).
The effect of Irinotecan on cancer is modeled by δ C Ir [ Ir ] [ C ] , because Irinotecan prevents DNA replication by inhibiting the topoisomerase 1 gene (TOP1) causing subsequent cell death [101,102,103]. Although used as a part of the FOLFIRI regimen, it is also given individually in the treatment of cancer.
Consequently, the equation for necrotic cells becomes:
d [ N ] d t = α N C ( δ C G β [ G β ] + δ C I γ [ I γ ] + δ C T C [ T C ] + δ C + δ C 5 fu [ 5 FU ] + δ C 5 fu I γ [ 5 FU ] [ I γ ] δ 5 fu M [ 5 FU ] [ M ] + δ C LV 5 fu [ 5 FU ] [ LV ] + δ C Ir [ Ir ] ) [ C ] δ N [ N ]

2.1.2. T-Cells

Within the tumor-microenvironment, 5-FU is known to help in the activation of T-cells, along with dendritic cells [104], and Irinotecan depletes the number of T-reg cells [105]. As reported in [104], the interaction between 5-FU, dendritic cells and T-cells is rather complex, but there is an overall increase in the activation and generation of helper T-cells and cytotoxic cells due to the indirect transfection of dendritic cells by 5-FU, modelled by introducing activation rates λ T h D 5 f u and λ T C D 5 f u into the Equations (2) and (3), respectively.
d [ T h ] d t = ( λ T h D [ D ] + λ T h D 5 fu [ 5 FU ] [ D ] + λ T h M [ M ] + λ T h μ 1 [ μ 1 ] ) [ T N ] ( δ T h μ 2 [ μ 2 ] + δ T h T r [ T r ] + δ T h ) [ T h ]
d [ T C ] d t = ( λ T C T h [ T h ] + λ T C D [ D ] + λ T C D 5 fu [ 5 FU ] [ D ] ) [ T N ] ( δ T C μ 2 [ μ 2 ] + δ T C T r [ T r ] + δ T C ) [ T C ]
In [48], activation rates for T-cells were made proportional to the density of naive T-cells in order to help stabilize the system. Therefore, we also add the above activation rates to the Equation (1).
d [ T N ] d t = A T N ( λ T h D [ D ] + λ T h D 5 fu [ 5 FU ] [ D ] + λ T h M [ M ] + λ T h μ 1 [ μ 1 ] ) [ T N ] ( λ T C T h [ T h ] + λ T C D [ D ] + λ T C D 5 fu [ 5 FU ] [ D ] ) [ T N ] ( λ T r T h [ T h ] + λ T r μ 2 [ μ 2 ] + λ T r G β [ G β ] ) [ T N ] δ T N [ T N ]
Among the major effects of Irinotecan on the tumor microenvironment, it is also shown to be the reduction of the abundance of regulatory T-cells [105]. There are other sources that report a significant reduction in T-regs after the chemotherapy using FOLFIRI [106].
d [ T r ] d t = ( λ T r T h [ T h ] + λ T r μ 2 [ μ 2 ] + λ T r G β [ G β ] ) [ T N ] ( δ T r μ 1 [ μ 1 ] + δ T r + δ T r Ir [ Ir ] ) [ T r ]

2.1.3. 5-FU & Leucovorin

Thus, combining the above individual interactions with 5-FU, the overall effect on the 5-FU concentration can be modeled by
d [ 5 F U ] d t = A i n j 5 f u ( t ) α 5 f u ( δ C 5 f u [ 5 F U ] + δ C 5 f u I γ [ 5 F U ] [ I γ ] δ 5 f u M [ 5 F U ] [ M ] + δ C L V 5 f u [ 5 F U ] [ L V ] ) [ C ] δ 5 f u D [ 5 F U ] [ D ] δ 5 f u [ 5 F U ]
Along with the parameters from the equation for cancer cells, 5-FU is modelled with additional parameters, δ 5 f u D for the amount of 5-FU used to help dendritic cells activate T helper cells, δ 5 f u to represent the elimination rate (about 80% of 5-FU is consumed in the liver [28,107,108]), α 5 f u as the amount of 5-FU used to kill a unit of cancer cells and A i n j 5 f u ( t ) is a function in time to model the repetitive cycle of 5FU dosages.
d [ L V ] d t = A i n j L V ( t ) δ L V [ L V ] α L V δ C L V 5 f u [ C ] [ 5 F U ] [ L V ]
Leucovorin is similarly modelled with A i n j L V ( t ) , a function for the dosage intake, a natural decay rate represented by δ L V , δ C L V 5 f u from the cancer equation since Leucovorin increases the cytotoxic effect of 5-FU and α L V the effectiveness of Leucovorin in killing cancer cells.

2.1.4. Irinotecan

Dynamics of Irinotecan is modeled in a similar way with a natural decay rate of Irinotecan denoted by δ I r , the death rate of cancer cells by Irinotecan as δ C I r and the effectiveness of killing cancer cells by the constant α C I r . Since we also know that Irinotecan depletes T-reg cells [105], we include the parameters δ T r I r , the depletion of T-reg cells by Irinotecan, and α I r T r , the effectiveness of Irinotecan in the depletion of T-reg cells.
d [ I r ] d t = A i n j I r ( t ) α I r C δ C I r [ C ] [ I r ] α I r T r δ T r I r [ I r ] [ T r ] δ I r [ I r ]

2.2. Non-Dimensionalization

Non-dimensionalization is used for additional numerical stability and to eliminate scale dependence [48]. The original system in [48] was non-dimensionalized by considering a non-dimensional variable X ¯ such that,
X ¯ = X X
for each variable X, where X is its steady-state value. For the new variables, namely the FOLFIRI drugs, we introduce new non-dimensional variables in the form of D ¯ for each drug D, defined as:
D ¯ = δ D [ D ] A i n j D d a i l y m e d i a n ,
where δ D is its natural decay rate and A i n j D d a i l y m e d i a n is its daily median dose per cycle from patient data (See Appendix A for further details).

2.3. Data of the Model

There are several popular tumor deconvolution methods to estimate immune cell frequencies using the gene expression profile of the tumors, and it has been shown in recent studies that CIBERSORTx method [109] has the highest accuracy among these methods [110,111,112]. In this study, we download RSEM normalized RNA-seq gene expression profiles in log 2 scale of the primary tumors of the 329 colon cancer patients from the TCGA project of COAD from UCSC Xena web portal [113] and transform it to the linear space. Then, we apply CIBERSORTx B-mode on the gene expression data to estimate immune cells frequencies.
In our previous study, K-means clustering of colon tumors based on their immune cells’ frequencies indicates that there are five distinct immune patterns of colonic tumors [48]. In this paper, we investigate the effect of FOLFIRI drugs on the same five distinct clusters. In each cluster, average immune cell frequencies that we use in the dynamical model have been shown in Figure 2, where the vertical bars show the standard deviations.
The data used in Figure 2 only gives us the ratio of immune cells so that we download TCGA biospecimen data from GDC portal that includes tumor dimension and necrotic cell percentage of the tumors to obtain the values of the model variables as described below.
Assuming the average amount of cancer cells is twice the average amount of total immune cells (TIC) and using the given necrosis percentage from TCGA biospecimen data, the average ratio of immune cells: cancer cells: necrotic cells is approximately 0.3:0.6:0.1 for colon cancer tumors. Also, we define size of tumor for each patient (P) as the product of the longest and the shortest dimension of the tumor, and total cell density (TCD) is assumed to be proportional to the size of the tumor.
TCD P = α d i m tumor size ( P ) 1 K all P tumor size ( P )
Using TCD and necrotic cell percentage ( N p ), we calculate the value of N and C in the following way:
N = TCD . N p , C = 2 3 TCD ( 1 N p ) and TIC = 0.5 C .
For scaling factor α d i m , we choose 7.5 × 10 4 to approximately match the average density of cancer cells for all patients to be 4.5 × 10 4 cells/cm 3 , which is reported in [114]. We determine the smallest tumors in each cluster and use their values as the initial conditions of the model for each cluster. As we solve the non-dimensionalized system, Table 2 shows the initial values of the non-dimensionalized variables, i.e., X / X for each cluster.

Treatment Data

For modeling FOLFIRI, we have downloaded the clinical drug data from the GDC portal that includes prescribed treatment information such as drug dosages, number of cycles of a specific treatment the patient received, days to start treatment and days to end treatment. To model each cluster, we have used patients average prescribed treatment information (Table 3). Note that minimum dosages values are used in parameter estimation of the model explained with more details in Appendix A.

2.4. Numerical Methods

To solve the system and study the dynamics of the variables, the previously developed code in [48] is modified to include the new equations, variables and parameters. The values for the parameters of the ODE system are either derived from research or by making appropriate assumptions based on biological information. The code uses SciPy solve_ivp function in python [115] to solve the system and the drug information is obtained from the treatment data as given in Table 3, to be used as the initial infusion rates. The infusion step function is modelled after the FOLFIRI chemotherapy regimen given in [116]. FOLFIRI is administered first with one hour drip of Irinotecan followed by one hour drip of Leucovorin, after which 5-FU is continuously infused into the bloodstream for 46 h. This treatment is repeated for the designated ‘number of cycles’, each cycle lasting over the period of the ‘cycle length’.

2.5. Sensitivity Analysis

To analyze the effect of parameter values on the dynamics of the system, we perform sensitivity analysis [117,118,119]. For the system d X d t = F X , θ , t consider (first order) sensitivity S of non-dimensional solution X with respect to the model parameters θ = θ i i = 1 , N ¯ to be defined as a vector
S i = d X d θ i , i = 1 , N ¯ .
As within the introduced model, the effects of the treatment do not extend to the steady state. Therefore, we consider time-dependent sensitivity satisfying the equation
d S i d t = F θ i + F X S i .
Additionally, we look at the “relative” sensitivity given by the formula
S ¯ i ( t ) = S i ( t ) θ i X ( t ) .
Relative sensitivity approach is commonly used in metabolic control analysis for biological reaction networks [120]. Then, for finite time T, we consider average sensitivity of each type:
S i = 1 T 0 T S i ( t ) d t , S ¯ i = 1 T 0 T S ¯ i ( t ) d t .

3. Results

3.1. Dynamics

We study the individual dynamics of each of the 17 variables and ‘Total cells’ (Immune cells, cancer and necrotic cells), with most of the parameter values obtained from [48] (using steady-state assumptions) and the new 19 parameters derived from parameter assumptions (as given in Appendix A). The drug inputs are the median drug dosages from patients treated with FOLFIRI—770, 725 and 300 mg of 5-FU, Leucovorin and Irinotecan, respectively. The average number of cycles and cycle length respectively are 12 and 14, as given in Table 3. In Figure 3, the initial conditions for the cells in each cluster are taken to be their steady-state values from [48], which correspond to the values of large tumors in each cluster.
Here, we study two different periods in the dynamics, during treatment and about two years after treatment has been stopped. We can observe that during the treatment, there is a decline in the number of Naive T-cells, T-reg cells, TGF- β and cancer cells. On the other hand, helper T-cells, cytotoxic cells, necrotic cells, HMGB1 and IFN- γ increase during treatment. Macrophages and μ 1 group of cytokines increase only by a small amount during treatment while μ 2 group of cytokines decrease slightly and this effect is most prominent for cluster 3 as compared to the other clusters. The effect of the treatment on both naive and active dendritic cells is not distinctly discernible, although we are able to observe smaller spikes during the time of treatment.
Based on our ODE system, we can infer that the increase in helper T-cells and cytotoxic cells is due to the activation rates λ T h 5 f u D and λ T C 5 f u D , that are modelled after the indirect activation of T-cells by dendritic cells and 5-FU [104]. On the other hand, the decline in T-reg cells must be due to the inhibitory effect of Irinotecan on T-reg cells, modelled in Equation (20) by the parameter δ T r I r . We have described macrophages to have an inhibitory effect on the cytotoxic action of 5-FU, but they themselves are not affected by this interaction, because this is the only interaction modelled between macrophages and the drugs in our network (as seen in Figure 1). Note, the small spikes of macrophages may be attributed to the activation rate of macrophages by T-helper cells, λ M T h , as 5-FU has an interaction with T-helper cells. Since we include the activation rates for helper and cytotoxic T-cells in the equation for naive T-cells (Equation (19)), the decrease in the density of naive cells is also reasonable. The small effect on dendritic cells, both naive and active can be explained by the parameter λ D C , since this parameter is linked to the density of cancer cells. For all the clusters, there is a rapid decline in the cancer cell density during treatment. Cluster 3 reaches the minimum cancer value among the clusters which is quite close to zero. Since the amount of necrotic cells increases with the decrease in cancer cells, it is reasonable to see it increases during treatment. Moreover, IFN- γ is not directly affected by any of the drugs, however, since helper T-cells increase during treatment this must also increase the value of IFN- γ . Cytokines, μ 1 , μ 2 , HMGB1 and TGF- β are all indirectly affected as per their interactions with macrophages, T-cells and dendritic cells in the tumor network Figure 1.
After the treatment stops, all the cells continue growing as per the model without any treatment, which is basically the system in [48], eventually reaching steady state. While the clusters generally show a common pattern, we are able to observe some unique behaviour in certain cell dynamics. Cluster 3 shows a slight increase in helper T-cells, cytotoxic cells, dendritic cells, macrophages, μ 1 , μ 2 and TGF- β , right after treatment and remains almost flat in its growth until the end of time period, decreasing again by a small amount. These combined changes could explain the decline in naive T-cells, as can be observed through Equation (19), and subsequently the decline in T-reg cells modeled by Equation (20), which depends on naive T-cells. Cluster 3 also reaches the lowest density of cancer cells, necrotic cells and ‘Total cells’ among all the clusters at the lowest point as well as by the end of the time period. While cluster 2 has the highest number of cancer cells and ‘Total cells’ at the lowest point and by the end of the time period. Importantly, comparing our results to that in [48], we can observe that cluster 3 shows a unique behavior for the model without treatment. In Figures 4 and 5 in [48], cluster 3 also shows similar patterns to that observed in our model with treatments, which further confirms that while the cell dynamics change during treatment, it does not change the overall dynamics.

3.2. Different Treatment Options

We investigate the individual and combined effect of the three drugs by plotting the cell dynamics with different combinations of the drugs. The initial conditions for the cells in each cluster are taken to be the steady state values, i.e., the values of large tumors [48]. For Figure 4A, we study the individual effect of 5-FU by keeping the other two drug values at zero. For Figure 4B, we study the combined effect of 5-FU and Leucovorin by keeping the Irinotecan dose at zero. For Figure 4C, we now study the combination of 5-FU and Irinotecan, keeping Leucovorin dose at zero. Finally, in Figure 4D, we study the combined effect of the three FOLFIRI drugs together. The drug values are their corresponding median dose values as given in Table 3. The number of cycles and cycle length respectively are 12 and 14 for all the drug combinations (given in Table 3).
There is not much differences in the dynamics of key players when using 5-FU alone versus 5-FU plus Leucovorin, except for a clear decrease only during and right after the treatments in cancer cells and consequently ‘Total cells’, which is to be expected since Leucovorin aids 5-FU in killing cancer cells but does not have any other direct interaction with other cells in the tumor microenvironment, as can be seen in Figure 1. With the addition of Irinotecan to 5-FU, there is a dramatic change in most of the cells, and especially T-reg cells. Adding Irinotecan changes the dynamics of T-reg and dendritic cells leading to slower tumor recurrence, especially for tumors in cluster 3. This is also an expected result, since Irinotecan is assumed to be 40% efficient in killing cancer cells, which is much higher than the other two drugs. Irinotecan also specifically targets T-reg cells, as is modelled by the parameter δ T r I r in the Equation (20).

3.3. Varying Treatment Start Time

We also investigate the effect of the treatment start time on dynamics of immune and cancer cells, since different patients start chemotherapy at different stages of their cancer. The initial conditions for the cells and cytokines for each cluster is obtained from taking the treatment data for the patient with the smallest tumor in each cluster (Table 2). The initial conditions for the drugs are taken to be their median dosages and the average for the number of cycles and cycle length as given in Table 3.
Giving the treatment at a later time, delays the time it takes for the cells to reach their steady states. As discussed previously, the cell dynamics are effected only during treatment, but the treatment does not change the overall dynamics, and we can observe the same effect here. In Figure 5A, cluster 5 has the lowest number of cancer cells after treatment, while in the other Figure 5B–D when the treatment is given after 3 years or later, cluster 3 reaches the lowest number of cancer cells after treatment. In Figure 5C,D, all cells seem to have already reached or very close to reaching their steady-states. Hence, there is not much of a difference between the dynamics in these cells between the plots C and D, except that the effect of the treatment is seen at a later time. In Figure 5D, since the treatment is given after about 7 years all the cells, including cancer and necrotic cells have already reached their steady states and therefore the dynamics is identical to that in Figure 3. As also observed in Figure 3 and Figure 4, cluster 3 shows some unique behavior compared to the other clusters. The number of dendritic cells in the other clusters is much higher than their initial value, but for cluster 3 it decreases after treatment and reaches to the same steady state value as cluster 5. Cluster 3 has the lowest number of total cells after the treatment in all plots, but cluster 5 has the lowest steady state value for cancer, necrotic and total cells.

3.4. Validating Model Using Patient Data

We have used clinical follow up data of the colon cancer patients downloaded from the GDC portal that includes tumor status and corresponding days of the follow up in different times to see if the number of cancer cells of the patients obtained from our model would match with their follow up data. Not every patients’ treatment information such as drug dosages and number of cycles are available, and there are only 4 patients who used Fluorouracil, Leucovorin and Irinotecan that we know their prescribed treatment information. For this reason, we validate our model with 36 patients that use Fluorouracil and Leucovorin. Note that most of these patients have also used other drugs such as Oxaliplatinum that are not included in the model.
We validate our model based on patients’ tumor status in their two different follow up days (Table 4). We consider patients who had only one follow up date in both groups. We exclude patients in the early follow up group, if their treatments have not been ended before their follow up day.

3.4.1. Predicting Tumor Status—ROC Curve

Data includes the start date of treatments and the last follow up date that the patients’ tumor status has been recorded. Therefore, to validate the model, we use the predicted numbers of cancer cells from the model at the exact number of days after each patient’s start date of treatments that their tumor status has been recorded. We then define a range of threshold values using minimum and maximum predicted values to create a ROC curve and predict patients’ tumor status as tumor free if the estimated number of values are less than the threshold value. As seen in Figure 6, area under the curve (AUC) for the first follow up days ROC curve (Figure 6A) is greater than the one for the last follow up days ROC curve (Figure 6B). As this Figure shows, the model’s predictions are good for the first follow up dates, but not for the second follow up dates. The bad prediction of the model for the last follow up date might be due to the fact that these patients might have had some other treatments, including surgery between the first and the last time of follow-up.

3.4.2. Individual Patients

We investigate the effect of dosages on the number of cancer and immune cells. From the same data set of 36 patients who were treated with 5-FU and Leucovorin, we choose one patient from each cluster. If the patient is ‘tumor-free’ on a particular follow-up day, then we expect that in our model results, the cancer on that day is at least less than the initial cancer value for small and medium size tumors and half of their initial values for large tumors. While for a ‘with-tumor’ patient, the number of cancer cells on the follow-up day should be greater than their initial values for small or medium size tumors or half of their initial values for large tumors. We find 10 patients satisfying these conditions among 15 patients from cluster 1, one patient among the 5 patients from each of the clusters 2 and 3 respectively, and 5 patients among the 10 patients from cluster 5. Note, there is no follow-up information available for cluster 4 and hence, cluster 4 patients have been excluded from the data. We investigate the effect of varying the dose on cancer and immune cells for some of patients that model predicts well to provide alternative optimal treatment options for these patients.
The drug dosages and initial conditions for the cells used in the model are based on the individual patient’s treatment data. In Figure 7, sub-figures A are plotted by keeping the 5-FU dose constant at the patient’s administered dose, while Leucovorin is varied between a range of 0.1 and 10 times the prescribed dose. Sub-figures B are plotted by keeping the Leucovorin dose at zero while varying the 5-FU dose between a range of 0.1 and 10 times the patient’s prescribed dose.
Among the patients in cluster 1, we investigate the results for the patient ‘TCGA-CM-6172’, since we know from the treatment data that this patient was treated with only 5-FU and Leucovorin without the aid of any other drug. Therefore, this patient is an ideal candidate for validating our model. This patient was administered a 662.5 mg dose of 5-FU and 574.17 mg dose of Leucovorin in 12 cycles and is reported to be ‘tumor-free’ at both the first and last follow-up days. As can be observed in the Figure 7(A1), the cancer at both follow-up days is less than the initial cancer value and therefore the follow-up information matches our model results.
We investigate the effect of varying the dose for the patient ‘TCGA-A6-6141’ in cluster 2 whose follow-up data matches our model results. This patient is treated with a 850 mg dose of 5-FU and 408.3 mg of Leucovorin in 12 cycles. This patient has also been reported ‘tumor-free’ at the follow-up day, and this is consistent with our results shown in Figure 7(A2) obtained for cancer.
We plot the dynamics for the patient ‘TCGA-A6-6142’ in cluster 3 whose tumor status matches our model results. This patient was administered 628.75 mg of 5-FU and 420.67 mg of Leucovorin in 12 cycles. This patient is reported to be ‘tumor-free’ at the first follow-up day, but is reported to be ‘with-tumor’ at the last follow-up day. This can also be confirmed through Figure 7(A3), where for the prescribed dose, the cancer is lower than the initial value at the first follow up day but grows back again surpassing the initial value, and are much higher on the last follow-up day. The increase in the total immune cells with increasing doses is much higher for this patient. Figure 7(A3,B3) also indicate a significant delay in reaching steady state for cancer when the dose increases. From these sub-figures, we can also recommend a higher dose of Leucovorin at 2103.33 mg or 5-FU alone at 3143.76 mg, or higher to achieve ‘tumor-free’ results on the last follow-up day.
We investigate the dynamics of one ‘with-tumor’ patient from cluster 5. Patient ‘TCGA-G4-6303’ was administered a dose of 660 mg of 5-FU and Leucovorin in 6 cycles. Figure 7(A4) shows that at the prescribed dose, the cancer is higher than the initial value. Note, the patient’s prescribed dosage does not even reduce the number of cancer cells much during treatment. 5-FU alone, does not have any significant effect in reducing the cancer for this patient, we instead recommend a dose of 3300 mg of Leucovorin along with the original 5-FU dose, or higher to achieve ‘tumor-free’ results for a long period of time.
Although the cancer initially decreases with treatment, it grows back again after a period of time according to the tumor recurrence rate modelled into the system. And therefore as one can observe for the above patients, the steady state values are always higher than the initial value, which means that the patient may once again have cancer years after receiving treatment. However, depending on the age of the patient, the tumor might not reach a visible size in the patients’ life time. In general, the best outcome is to achieve a steady-state value for cancer that remains below the initial value. In order to demonstrate this, we plot the patient ‘TCGA-A6-6781’ from cluster 5, for whom, according our model, its steady state value is less than its initial value. This patient is treated with 902.67 mg of 5-FU and 166.67 mg of Leucovorin in 12 cycles. Although there is nothing distinctive to predict such results, in fact this patient has the highest initial cancer value among all the patients in the dataset. Interestingly, it seems 5-FU alone in a higher dose is a better treatment option for this patient, while for most patients the combination of 5-FU and Leucovorin works better than 5-FU alone.

3.4.3. Effect of Sensitive Parameters

We perform sensitivity analysis of the non-dimensional system with respect to treatment-related parameters. Cells are assumed to be at their “no treatment” steady states (i.e., large tumors) for the initial conditions. The resulting time-averaged sensitivity for most sensitive parameters is presented on Figure 8.
These results show that the most sensitive parameter is δ C 5 f u L V , and δ 5 f u M is the most sensitive among immune parameters. We investigate their effect on cancer cells and ‘Total cells’, by varying them between a range of 0.1 and 5 times their original value obtained from parameters’ assumptions, and plotting them against the number of cells at different time points. In these plots, the initial conditions for the cells were chosen to be their steady state values, i.e., values of large tumors.
Figure 9 shows a decline in the number of cancer cells with an increase in the δ C 5 f u L V value, and this decline is observed to be smaller in cluster 2 as compared with the other clusters. The lowest number of cancer cells can be observed at 169 days, which is right after the treatment ends. Total cells seem to only be effected by increasing values at the 2nd and 3rd year time points, while for the other time-points there is barely any effect. This must be due to the fact that after treatment ends, the cancer starts to grow back again. In clusters 1 and 5, although the total number of cells first decreases at the year three, it increases at higher δ C 5 f u L V values.
Varying δ 5 f u M does not seem to have any significant change in the cancer or total cells (Figure 10), which is consistent with noticeably lower sensitivity level of this parameter compared to δ C 5 f u L V (Figure 8). We similarly see the lowest number of cancer cells at 169 days, right after stopping the treatment, and the cancer increases after treatment.

4. Discussion

Cancer is a heterogeneous disease that includes different components such as immune cells, cancer cells or lymphatic vessels [121]. Drug resistance is one of the main problems in cancer studies [122]. Mechanisms or components of cancer are usually investigated one by one in traditional in vitro and in vivo studies. Although these studies provide valuable information about a mechanism, each of these studies alone is not able to provide necessary and sufficient information to explain cancer complexity [123]. Having more accessible biological experimental data set and new advances in tumor deconvolution methods lead to increase demand in data driven mathematical models that help us to model several mechanisms together and study the complexity of the system in a more effective way [124]. Tumor microenvironment components have an essential role in the explanation of poor prognosis and immune escape in CRC [21,22,23,125], and they have been used to explain chemotherapy and immunotherapy sensitivity in many studies [44,126,127]. In this study, we develop a data driven mathematical model that takes each tumor’s characteristic into consideration for the treatment. We have used interaction between TME components and drug mechanism to model the response to the FOLFIRI treatment in colon cancer using gene expression profiles of the CRC primary tumors to estimate immune patterns. Our results demonstrate how leucovorin increases the efficiency of 5-FU on cancer cells (Figure 4A,B) that has been shown in many studies [94,98,99]. Also, the impact of the combination of irinotecan, 5-FU, and leucovorin can be seen clearly on cancer cells in our model (Figure 4D) as reported in other studies [31,128]. In addition, we have applied treatment in different time points to show how cancer comes back aggressively if the treatment starts late (Figure 5).
The mathematical model demonstrates the relation between immune infiltration and drug’s effects on CRC primary tumors. Relative change in T-reg/T-helper ratio has been found as clinical index for response prediction; colon cancer patients with higher T-reg/Th ratio respond better to treatments [129]. Cluster 2 that has the lowest level of T-reg/T-helper ratio (Figure 2) is more resistant to FOLFIRI treatment since after the treatment the number of cancer cells increase faster and the values are not close to zero right after the treatment, while in other clusters the number of cancer cells approaches zero and their growth rate is slower (Figure 3). In contrast, the FOLFIRI treatment works better for cluster 3 that has the highest level of T-reg/T-helper ratio compared to the other clusters. Also, a decrease in regulatory T-cells after FOLFIRI treatment in colon cancer compared to pre-treatment level has been found to be associated with better survival months [130]. We observe a similar decrease for regulatory T-cells in tumors in cluster 3 that have a good response to the FOLFIRI treatment (Figure 3).
It has been observed that the number of T-reg cells significantly decreases for colon cancer patients who have a high number of regulatory T-cells before FOLFIRI treatment [106]. We divide patients into two groups based on their regulatory T-cells values as high T-reg and low T-reg. When we compare T-reg values before and after the treatment, we also see a decrease but it is not significant. In our validation data, there are only a few number of patients in the high T-reg group, and that might be the reason for not observing a significant p-value. It is important to mention that it would be ideal to use the gene expression of the patients after the treatment or data from patients who do not use other drugs rather than FOLFIRI to validate the model. However, follow up gene expression values are mostly unavailable, and we have only a small number of patients with all their treatment information available. Thus, we validate our model based on patients tumor status at follow up date. As we see in Figure 6, our model predicts much better for the first follow up data compare to the last follow up. That might be reasonable, because patients might have had other treatments such as surgery between their first and last time of follow ups.
In general, this work has some limitations that should be considered when these results are used. As mentioned above, patients included in the validation were also treated with other drugs such as oxaliplatin. Moreover, only 31 patients have early follow up data, and only 36 patients have late follow up data. Therefore, the validations have been done on a small number of patients who might have undergone other treatments. We have not had any data to validate the predictions of the model for the responses to Irinotecan.
Although our model has some limitation due to the lack of time course data, it presents valuable insight about the interactions between FOLFIRI treatment and the tumor micro-environment. Moreover, many studies can build upon on this one to provide the best treatment options for patients using only patients’ gene expression data. One way to improve this model is adding other chemotherapy options such as Oxaliplatinum, and different parameter fitting algorithm can be applied to increase the accuracy of the model [131,132,133,134]. Another possible way might be including other drug resistance factors in the model and extend it to a powerful model that considers other mechanisms and other types of patients’ data.

Author Contributions

Conceptualization, L.S.; methodology, A.B., S.S., A.K., L.S.; software, A.B., S.S., A.K.; validation, A.B., S.S.; formal analysis, A.B., S.S., A.K.; investigation, A.B., S.S., A.K.; resources, L.S.; data curation, S.S.; writing—original draft preparation, A.B., S.S., A.K., L.S.; writing—review and editing, L.S., S.S., A.B.; visualization, A.B., S.S.; supervision, L.S.; project administration, L.S.; funding acquisition, L.S. All authors have read and agreed to the published version of the manuscript.

Funding

Research reported in this publication was supported by the National Cancer Institute of the National Institutes of Health under Award Number R21CA242933. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The TCGA data of COAD project underlying this article are available at https://portal.gdc.cancer.gov and https://xenabrowser.net/datapages/ (accessed on 30 March 2021) [113]. Codes are available at our GitHub page https://github.com/ShahriyariLab/Data-driven-mathematical-model-of-FOLFIRI-treatment-for-colon-cancer (accessed on 30 March 2021).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
COADColon adenocarcinoma
CRCColorectal cancer
HMGB1High mobility group box 1
DAMPdamage-associated molecular pattern
TAMTumor-associated macrophages
ODEOrdinary differential equation
TYMSthymidylate synthase
TCDTotal cell density
TICTotal immune cell

Appendix A. ODE System & Non-Dimensionalization

The entire ODE system, which has been modeled in this paper is given below. The activation and proliferation rates are denoted by λ , while death and degradation rates are denoted by δ . Changes and additions to the system modeled in [48] are highlighted in bold.
d [ T N ] d t = A T N ( λ T h D [ D ] + λ T h D 5 fu [ 5 FU ] [ D ] + λ T h M [ M ] + λ T h μ 1 [ μ 1 ] ) [ T N ] ( λ T C T h [ T h ] + λ T C D [ D ] + λ T C D 5 fu [ 5 FU ] [ D ] ) [ T N ] ( λ T r T h [ T h ] + λ T r μ 2 [ μ 2 ] + λ T r G β [ G β ] ) [ T N ] δ T N [ T N ]
d [ T h ] d t = ( λ T h D [ D ] + λ T h D 5 fu [ 5 FU ] [ D ] + λ T h M [ M ] + λ T h μ 1 [ μ 1 ] ) [ T N ] ( δ T h μ 2 [ μ 2 ] + δ T h T r [ T r ] + δ T h ) [ T h ]
d [ T C ] d t = ( λ T C T h [ T h ] + λ T C D [ D ] + λ T C D 5 fu [ 5 FU ] [ D ] ) [ T N ] ( δ T C μ 2 [ μ 2 ] + δ T C T r [ T r ] + δ T C ) [ T C ]
d [ T r ] d t = ( λ T r T h [ T h ] + λ T r μ 2 [ μ 2 ] + λ T r G β [ G β ] ) [ T N ] ( δ T r μ 1 [ μ 1 ] + δ T r + δ T r Ir [ Ir ] ) [ T r ]
d D N d t = A D N λ D H H + λ D C C D N δ D H H + δ D D N
d D d t = λ D H H + λ D C C D N δ D H H + δ D C C + δ D D
d M d t = λ M μ 2 μ 2 + λ M I γ I γ + λ M T h T h M 0 M δ M M
d [ C ] d t = ( λ C + λ C μ 1 [ μ 1 ] ) [ C ] ( 1 [ C ] C 0 ) ( δ C G β [ G β ] + δ C I γ [ I γ ] + δ C T C [ T C ] + δ C + δ C 5 fu [ 5 FU ] + δ C 5 fu I γ [ 5 FU ] [ I γ ] δ 5 fu M [ 5 FU ] [ M ] + δ C LV 5 fu [ 5 FU ] [ LV ] + δ C Ir [ Ir ] ) [ C ]
d [ N ] d t = α N C ( δ C G β [ G β ] + δ C I γ [ I γ ] + δ C T C [ T C ] + δ C + δ C 5 fu [ 5 FU ] + δ C 5 fu I γ [ 5 FU ] [ I γ ] δ 5 fu M [ 5 FU ] [ M ] + δ C LV 5 fu [ 5 FU ] [ LV ] + δ C Ir [ Ir ] ) [ C ] δ N [ N ]
d H d t = λ H N N + λ H M M + λ H T h T h + λ H T C T C + λ H T r T r δ H H
d μ 1 d t = λ μ 1 T h T h + λ μ 1 M M + λ μ 1 D D δ μ 1 μ 1
d μ 2 d t = λ μ 2 M M + λ μ 2 D D + λ μ 2 T r T r δ μ 2 μ 2
d I γ d t = λ I γ T h T h + λ I γ T C T C + λ I γ M M δ I γ I γ
d G β d t = λ G β M M + λ G β T r T r δ G β G β
d [ 5 FU ] dt = A inj 5 fu ( t ) α 5 fu ( δ C 5 fu [ 5 FU ] + δ C 5 fu I γ [ 5 FU ] [ I γ ] δ 5 fu M [ 5 FU ] [ M ] + δ C LV 5 fu [ 5 FU ] [ LV ] ) [ C ] δ 5 fu D [ 5 FU ] [ D ] δ 5 fu [ 5 FU ]
d [ Ir ] dt = A inj Ir ( t ) α Ir C δ C Ir [ C ] [ Ir ] α Ir T r δ T r Ir [ Ir ] [ T r ] δ Ir [ Ir ]
d [ LV ] dt = A inj LV ( t ) δ LV [ LV ] α LV δ C LV 5 fu [ C ] [ 5 FU ] [ LV ]
We introduce non-dimensional variables,
5 F U ¯ = δ 5 f u [ 5 F U ] A i n j 5 f u d a i l y m e d i a n I r ¯ = δ I r [ I r ] A i n j I r d a i l y m e d i a n L V ¯ = δ L V [ L V ] A i n j L V d a i l y m e d i a n ,
then we have the following non-dimensional equations (only the new and altered ones have been included, the rest of non-dimensionalized equations are exactly the same as the ones provided in [48]):
d [ 5 FU ¯ ] dt = A ¯ inj 5 fu ( t ) α ¯ 5 fu ( δ ¯ C 5 fu [ 5 FU ¯ ] + δ ¯ C 5 fu I γ [ 5 FU ¯ ] [ I γ ¯ ] δ ¯ 5 fu M [ 5 FU ¯ ] [ M ¯ ] + δ ¯ C LV 5 fu [ 5 FU ¯ ] [ LV ¯ ] ) [ C ¯ ] δ ¯ 5 fu D [ 5 FU ¯ ] [ D ¯ ] δ 5 fu [ 5 FU ¯ ]
d [ Ir ¯ ] dt = A ¯ inj Ir ( t ) α ¯ Ir C δ ¯ C Ir [ C ¯ ] [ Ir ¯ ] α ¯ Ir T r δ ¯ T r Ir [ Ir ¯ ] [ T ¯ r ] δ Ir [ Ir ¯ ]
d [ LV ¯ ] dt = A ¯ inj LV ( t ) δ LV [ LV ¯ ] α ¯ LV δ ¯ C LV 5 fu [ C ¯ ] [ 5 FU ¯ ] [ LV ¯ ]
d [ T r ¯ ] d t = ( λ ¯ T r T h [ T ¯ h ] + λ ¯ T r μ 2 [ μ 2 ] + λ ¯ T r G β [ G ¯ β ] ) [ T ¯ N ] ( δ ¯ T r μ 1 [ μ ¯ 1 ] + δ T r + δ ¯ T r Ir [ Ir ¯ ] ) [ T ¯ r ]
d [ T h ¯ ] d t = ( λ ¯ T h D [ D ¯ ] + λ ¯ T h D 5 fu [ 5 FU ¯ ] [ D ¯ ] + λ ¯ T h M [ M ¯ ] + λ ¯ T h μ 1 [ μ ¯ 1 ] ) [ T N ¯ ] ( δ ¯ T h μ 2 [ μ ¯ 2 ] + δ ¯ T h T r [ T ¯ r ] + δ T h ) [ T h ¯ ]
d [ T C ¯ ] d t = ( λ ¯ T C T h [ T h ¯ ] + λ ¯ T C D [ D ¯ ] + λ ¯ T C D 5 fu [ 5 FU ¯ ] [ D ¯ ] ) [ T N ¯ ] ( δ ¯ T C μ 2 [ μ ¯ 2 ] + δ ¯ T C T r [ T ¯ r ] + δ T C ) [ T C ¯ ]
d [ T N ¯ ] d t = A ¯ T N α T N T h ( λ ¯ T h D [ D ¯ ] + λ ¯ T h D 5 fu [ 5 FU ¯ ] [ D ¯ ] + λ ¯ T h M [ M ¯ ] + λ ¯ T h μ 1 [ μ ¯ 1 ] ) [ T N ¯ ] α T N T C ( λ ¯ T C T h [ T h ¯ ] + λ ¯ T C D [ D ¯ ] + λ ¯ T C D 5 fu [ 5 FU ¯ ] [ D ¯ ] ) [ T N ¯ ] α T N T r ( λ ¯ T r T h [ T h ¯ ] + λ ¯ T r μ 2 [ μ ¯ 2 ] + λ ¯ T r G β [ G ¯ β ] ) [ T N ¯ ] δ T N [ T N ¯ ]
d [ C ] d t = ( λ C + λ ¯ C μ 1 [ μ ¯ 1 ] ) [ C ¯ ] 1 [ C ] C 0 ( δ ¯ C G β [ G ¯ β ] + δ ¯ C I γ [ I ¯ γ ] + δ ¯ C T C [ T ¯ C ] + δ C + δ ¯ C 5 fu [ 5 FU ¯ ] + δ ¯ C 5 fu I γ [ 5 FU ¯ ] [ I ¯ γ ] δ ¯ 5 fu M [ 5 FU ¯ ] [ M ¯ ] + δ ¯ C LV 5 fu [ 5 FU ¯ ] [ LV ¯ ] + δ ¯ C Ir [ Ir ¯ ] ) [ C ¯ ]
where the newly introduced non-dimensional parameters are:
A ¯ i n j 5 f u = A i n j 5 f u δ 5 f u A i n j 5 f u d a i l y m e d i a n α ¯ 5 f u = α 5 f u C δ 5 f u A i n j 5 f u d a i l y m e d i a n λ ¯ T h D 5 f u = λ T h D 5 f u D T N A i n j 5 f u d a i l y m e d i a n δ 5 f u T h δ ¯ 5 f u D = δ 5 f u D D λ ¯ T C D 5 f u = λ T C D 5 f u D T N A i n j 5 f u d a i l y m e d i a n δ 5 f u T C δ ¯ C 5 f u = δ C 5 f u A i n j 5 f u d a i l y m e d i a n δ 5 f u δ ¯ C 5 f u I γ = δ C 5 f u I γ I γ A i n j 5 f u d a i l y m e d i a n δ 5 f u δ ¯ 5 f u M = δ 5 f u M M A i n j 5 f u d a i l y m e d i a n δ 5 f u δ ¯ C L V 5 f u = δ C L V 5 f u · [ 5 F U ¯ ] A i n j 5 f u d a i l y m e d i a n δ 5 f u · [ L V ¯ ] A i n j L V d a i l y m e d i a n δ L V A ¯ i n j L V = A i n j L V δ L V A i n j L V d a i l y m e d i a n α ¯ L V = α L V C δ L V A i n j L V d a i l y m e d i a n δ ¯ C I r = δ C I r A i n j I r d a i l y m e d i a n δ I r A ¯ i n j I r = A i n j I r δ I r A i n j I r d a i l y m e d i a n δ ¯ T r I r = δ T r I r A i n j I r d a i l y m e d i a n δ I r α ¯ I r C = α I r C C δ I r A i n j I r d a i l y m e d i a n α ¯ I r T r = α I r T r T r δ I r A i n j I r d a i l y m e d i a n
Note that A i n j ( t ) represents the infusion step function and A i n j d a i l y m i n / m e d i a n is the corresponding drug minimum or median administered dose based on patients’ data as given in Table 3.
For the natural decay rates of the drugs, we use the same formula from [48] by using their respective terminal/elimination half-lives as reported in available studies.
From available studies, we know that the terminal half-life of 5-FU is in between 8 and 20 min [135]. Therefore, using the average half-life as 14 min, we have
δ 5 f u = ln 2 t 1 / 2 = ln 2 14 mins = 71.3 day 1 .
Since Leucovorin is a mixture of the two isomers [6R]-5-formyl-tetrahydrofolate and [6S]-5-formyl-tetrahydrofolate, both with very different half-lives, to determine its half-life, we use [136], which reports it as 5.7 h and the terminal half-life as given by [137] for the total folates and metabolites of 7.59 h. Therefore, taking the average half-life as 6.5 h (approx. 5–8 h) we have
δ L V = ln 2 6.5 hrs = 2.56 day 1
And the reported terminal half-life of Irinotecan ranges between 5 and 18 h [103,138,139], giving an average half-life of 11.5 h. Therefore,
δ I r = ln 2 11.5 hrs = 1.45 day 1
The median dosages are used for the parameters A i n j d a i l y m e d i a n corresponding to each drug by dividing them by the average cycle length to obtain the daily dose:
A i n j d a i l y m e d i a n 5 f u = m e d i a n d o s a g e p e r c y c l e p e r d a y = 770 14
A i n j d a i l y m e d i a n I r = m e d i a n d o s a g e p e r c y c l e p e r d a y = 300 14
A i n j d a i l y m e d i a n L V = m e d i a n d o s a g e p e r c y c l e p e r d a y = 725 14
Similarly, the minimum daily dose per cycle from patients’ data for each corresponding drug is the minimum dose divided by the average cycle length:
A i n j d a i l y m i n 5 f u = m i n . d o s a g e p e r c y c l e p e r d a y = 598 14
A i n j d a i l y m i n I r = m i n . d o s a g e p e r c y c l e p e r d a y = 208 14
A i n j d a i l y m i n L V = m i n . d o s a g e p e r c y c l e p e r d a y = 75 14
We used the following assumptions for the parameter values (all in dimensional form). Values for α s have been derived from their efficiency, either from available studies or by appropriate assumptions. Note that maximum values and steady state values of the variables are taken from [48].
δ 5 f u I γ · I γ m a x I γ = δ C 5 f u · 3 2 δ 5 f u M · M m a x M = δ ¯ 5 f u D · D m a x D = 1 2 · δ C A i n j 5 f u d a i l y m i n δ C 5 f u · A i n j 5 f u d a i l y m i n = 10 2 · δ C λ T h D 5 f u · A i n j 5 f u d a i l y m i n = 1 2 · λ ¯ T h D λ T C D 5 f u · A i n j 5 f u d a i l y m i n = 1 2 · λ ¯ T C D α 5 f u = derived from efficiency
It is assumed that about 20% of 5-FU is effectively used to kill cancer cells (as reported above [28,107,108]).
δ C 5 f u L V = 0.1 · δ C 5 f u A i n j L V d a i l y m i n α L V = derived from efficiency
assuming 5% efficiency of killing cancer cells for Leucovorin.
δ C I r = 0.5 · δ C 5 f u · A i n j 5 f u d a i l y m i n A i n j I r d a i l y m i n
δ I r T r = δ C I r
α I r T r = α I r C
α I r C = derived from efficiency
assuming 40% efficiency of killing cancer cells for Irinotecan, since at least 60% of the drug is known to be eliminated from the body without being used [103].
In non-dimensional form, the parameter assumptions are as follows:
λ ¯ T h D 5 f u = λ ¯ T h D A i n j 5 f u d a i l y m e d i a n 2 · δ 5 f u A i n j 5 f u d a i l y m i n λ ¯ T C D 5 f u = λ ¯ T C D A i n j 5 f u d a i l y m e d i a n 2 · δ 5 f u A i n j 5 f u d a i l y m i n δ ¯ C 5 f u I γ = 3 · δ ¯ C 5 f u · I γ 2 · I γ m a x δ ¯ 5 f u M = δ C · M · A i n j 5 f u d a i l y m e d i a n 2 · δ 5 f u · M m a x · A i n j 5 f u d a i l y m i n δ ¯ 5 f u D = δ C · D · A i n j 5 f u d a i l y m e d i a n 2 · δ 5 f u · D m a x · A i n j 5 f u d a i l y m i n δ ¯ C 5 f u = 10 2 · A i n j 5 f u d a i l y m e d i a n A i n j 5 f u d a i l y m i n · δ C δ 5 f u δ ¯ C I r = 1 2 · A i n j I R d a i l y m e d i a n A i n j I R d a i l y m i n δ 5 f u δ ¯ C 5 f u δ I R A i n j I R d a i l y m i n A i n j I R d a i l y m e d i a n δ ¯ T r I r = δ ¯ C I r δ ¯ C 5 f u L V = 0.1 · A i n j 5 f u d a i l y m e d i a n A i n j L V d a i l y m i n · δ ¯ C 5 f u δ L V α ¯ L V = L V e f f 1 L V e f f · δ L V δ ¯ C 5 f u L V ( 1 5 f u e f f ) α ¯ I r C = I R e f f 1 I R e f f δ I r δ ¯ C I r α ¯ I r T r = I R e f f 1 I R e f f δ I r δ ¯ T r I r α ¯ 5 f u = 5 F U e f f 1 5 F U e f f δ 5 f u δ ¯ 5 f u D δ ¯ C 5 f u + δ ¯ C 5 f u I γ δ 5 f u M + δ C 5 f u L V · ( 1 L V e f f )

References

  1. Key Statistics for Colorectal Cancer. Available online: https://www.cancer.org/cancer/colon-rectal-cancer/about/key-statistics.html (accessed on 30 March 2021).
  2. Guinney, J.; Dienstmann, R.; Wang, X.; de Reyniès, A.; Schlicker, A.; Soneson, C.; Marisa, L.; Roepman, P.; Nyamundanda, G.; Angelino, P.; et al. The consensus molecular subtypes of colorectal cancer. Nat. Med. 2015, 21, 1350–1356. [Google Scholar] [CrossRef]
  3. Van Cutsem, E.; Cervantes, A.; Adam, R.; Sobrero, A.; Van Krieken, J.; Aderka, D.; Aranda Aguilar, E.; Bardelli, A.; Benson, A.; Bodoky, G.; et al. ESMO consensus guidelines for the management of patients with metastatic colorectal cancer. Ann. Oncol. 2016, 27, 1386–1422. [Google Scholar] [CrossRef]
  4. Poston, G.J.; Tait, D.; O’Connell, S.; Bennett, A.; Berendse, S. Diagnosis and management of colorectal cancer: Summary of NICE guidance. BMJ 2011, 343, d6751. [Google Scholar] [CrossRef] [PubMed]
  5. Siegel, R.L.; Miller, K.D.; Jemal, A. Cancer statistics, 2020. CA Cancer J. Clin. 2020, 70, 7–30. [Google Scholar] [CrossRef] [PubMed]
  6. Morse, M.A.; Hochster, H.; Benson, A. Perspectives on Treatment of Metastatic Colorectal Cancer with Immune Checkpoint Inhibitor Therapy. Oncologist 2020, 25, 33–45. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Wu, T.; Dai, Y. Tumor microenvironment and therapeutic response. Cancer Lett. 2017, 387, 61–68. [Google Scholar] [CrossRef]
  8. Golden, E.B.; Frances, D.; Pellicciotta, I.; Demaria, S.; Helen Barcellos-Hoff, M.; Formenti, S.C. Radiation fosters dose-dependent and chemotherapy-induced immunogenic cell death. OncoImmunology 2014, 3, e28518. [Google Scholar] [CrossRef] [Green Version]
  9. Schildkopf, P.; Frey, B.; Mantel, F.; Ott, O.J.; Weiss, E.M.; Sieber, R.; Janko, C.; Sauer, R.; Fietkau, R.; Gaipl, U.S. Application of hyperthermia in addition to ionizing irradiation fosters necrotic cell death and HMGB1 release of colorectal tumor cells. Biochem. Biophys. Res. Commun. 2010, 391, 1014–1020. [Google Scholar] [CrossRef] [PubMed]
  10. Apetoh, L.; Ghiringhelli, F.; Tesniere, A.; Criollo, A.; Ortiz, C.; Lidereau, R.; Mariette, C.; Chaput, N.; Mira, J.P.P.; Delaloge, S.; et al. The interaction between HMGB1 and TLR4 dictates the outcome of anticancer chemotherapy and radiotherapy. Immunol. Rev. 2007, 220, 47–59. [Google Scholar] [CrossRef]
  11. Liu, L.; Yang, M.; Kang, R.; Wang, Z.; Zhao, Y.; Yu, Y.; Xie, M.; Yin, X.; Livesey, K.M.; Lotze, M.T.; et al. HMGB1-induced autophagy promotes chemotherapy resistance in leukemia cells. Leuk. Off. J. Leuk. Soc. Am. Leuk. Res. Fund UK 2011, 25, 23–31. [Google Scholar] [CrossRef] [Green Version]
  12. Scaffidi, P.; Misteli, T.; Bianchi, M.E. Release of chromatin protein HMGB1 by necrotic cells triggers inflammation. Nature 2002, 418, 191–195. [Google Scholar] [CrossRef]
  13. Lotze, M.T.; Tracey, K.J. High-mobility group box 1 protein (HMGB1): Nuclear weapon in the immune arsenal. Nat. Rev. Immunol. 2005, 5, 331–342. [Google Scholar] [CrossRef]
  14. Xu, X.; Fu, X.Y.; Plate, J.; Chong, A.S. IFN-gamma induces cell growth inhibition by Fas-mediated apoptosis: Requirement of STAT1 protein for up-regulation of Fas and FasL expression. Cancer Res. 1998, 58, 2832–2837. [Google Scholar] [PubMed]
  15. Kroemer, G.; Galluzzi, L.; Kepp, O.; Zitvogel, L. Immunogenic Cell Death in Cancer Therapy. Annu. Rev. Immunol. 2013, 31, 51–72. [Google Scholar] [CrossRef] [PubMed]
  16. Nathan, C.F.; Murray, H.W.; Wiebe, M.E.; Rubin, B.Y. Identification of interferon-gamma as the lymphokine that activates human macrophage oxidative metabolism and antimicrobial activity. J. Exp. Med. 1983, 158, 670–689. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Bogdan, C.; Stenger, S.; Röllinghoff, M.; Solbach, W. Cytokine interactions in experimental cutaneous leishmaniasis. Interleukin 4 synergizes with interferon-γ to activate murine macrophages for killing ofLeishmania major amastigotes. Eur. J. Immunol. 1991, 21, 327–333. [Google Scholar] [CrossRef] [PubMed]
  18. Popivanova, B.K.; Kitamura, K.; Wu, Y.; Kondo, T.; Kagaya, T.; Kaneko, S.; Oshima, M.; Fujii, C.; Mukaida, N. Blocking TNF-alpha in mice reduces colorectal carcinogenesis associated with chronic colitis. J. Clin. Investig. 2008, 118, 560–570. [Google Scholar] [CrossRef] [PubMed]
  19. Qian, B.Z.; Pollard, J.W. Macrophage Diversity Enhances Tumor Progression and Metastasis. Cell 2010, 141, 39–51. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Ostuni, R.; Kratochvill, F.; Murray, P.J.; Natoli, G. Macrophages and cancer: From mechanisms to therapeutic implications. Trends Immunol. 2015, 36, 229–239. [Google Scholar] [CrossRef]
  21. Pagès, F.; Berger, A.; Camus, M.; Sanchez-Cabo, F.; Costes, A.; Molidor, R.; Mlecnik, B.; Kirilovsky, A.; Nilsson, M.; Damotte, D.; et al. Effector memory T cells, early metastasis, and survival in colorectal cancer. N. Engl. J. Med. 2005, 353, 2654–2666. [Google Scholar] [CrossRef]
  22. Dahlin, A.M.; Henriksson, M.L.; Van Guelpen, B.; Stenling, R.; Öberg, Å.; Rutegård, J.; Palmqvist, R. Colorectal cancer prognosis depends on T-cell infiltration and molecular characteristics of the tumor. Mod. Pathol. 2011, 24, 671–682. [Google Scholar] [CrossRef]
  23. Funada, Y.; Noguchi, T.; Kikuchi, R.; Takeno, S.; Uchida, Y.; Gabbert, H.E. Prognostic significance of CD8+ T cell and macrophage peritumoral infiltration in colorectal cancer. Oncol. Rep. 2003, 10, 309–313. [Google Scholar] [CrossRef]
  24. Lugade, A.A.; Sorensen, E.W.; Gerber, S.A.; Moran, J.P.; Frelinger, J.G.; Lord, E.M. Radiation-Induced IFN-γ Production within the Tumor Microenvironment Influences Antitumor Immunity. J. Immunol. 2008, 180, 3132–3139. [Google Scholar] [CrossRef]
  25. Burnette, B.C.; Liang, H.; Lee, Y.; Chlewicki, L.; Khodarev, N.N.; Weichselbaum, R.R.; Fu, Y.X.; Auh, S.L. The Efficacy of Radiotherapy Relies upon Induction of Type I Interferon–Dependent Innate and Adaptive Immunity. Cancer Res. 2011, 71, 2488–2496. [Google Scholar] [CrossRef] [Green Version]
  26. Tosolini, M.; Kirilovsky, A.; Mlecnik, B.; Fredriksen, T.; Mauger, S.; Bindea, G.; Berger, A.; Bruneval, P.; Fridman, W.H.; Pagès, F.; et al. Clinical impact of different classes of infiltrating T cytotoxic and helper cells (Th1, th2, treg, th17) in patients with colorectal cancer. Cancer Res. 2011, 71, 1263–1271. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Mantovani, A.; Allavena, P. The interaction of anticancer therapies with tumor-associated macrophages. J. Exp. Med. 2015, 212, 435–445. [Google Scholar] [CrossRef] [PubMed]
  28. Thorn, C.F.; Marsh, S.; Carrillo, M.W.; McLeod, H.L.; Klein, T.E.; Altman, R.B. PharmGKB summary: Fluoropyrimidine pathways. Pharmacogenet. Genom. 2010, 21, 237–242. [Google Scholar] [CrossRef] [Green Version]
  29. Longley, D.B.; Harkin, D.P.; Johnston, P.G. 5-Fluorouracil: Mechanisms of action and clinical strategies. Nat. Rev. Cancer 2003, 3, 330–338. [Google Scholar] [CrossRef]
  30. Johnston, P.G.; Kaye, S. Capecitabine: A novel agent for the treatment of solid tumors. Anti-Cancer Drugs 2001, 12, 639–646. [Google Scholar] [CrossRef] [PubMed]
  31. Douillard, J.; Cunningham, D.; Roth, A.; Navarro, M.; James, R.; Karasek, P.; Jandik, P.; Iveson, T.; Carmichael, J.; Alakl, M.; et al. Irinotecan combined with fluorouracil compared with fluorouracil alone as first-line treatment for metastatic colorectal cancer: A multicentre randomised trial. Lancet 2000, 355, 1041–1047. [Google Scholar] [CrossRef]
  32. Shahriyari, L.; Komarova, N.L. Symmetric vs. asymmetric stem cell divisions: An adaptation against cancer? PLoS ONE 2013, 8, e76195. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Shahriyari, L.; Komarova, N.L. The role of the bi-compartmental stem cell niche in delaying cancer. Phys. Biol. 2015, 12, 055001. [Google Scholar] [CrossRef]
  34. Shahriyari, L.; Mahdipour-Shirayeh, A. Modeling dynamics of mutants in heterogeneous stem cell niche. Phys. Biol. 2017, 14. [Google Scholar] [CrossRef]
  35. Bollas, A.; Shahriyari, L. The role of backward cell migration in two-hit mutants’ production in the stem cell niche. PLoS ONE 2017, 12, e0184651. [Google Scholar] [CrossRef] [Green Version]
  36. Brady, R.; Enderling, H. Mathematical Models of Cancer: When to Predict Novel Therapies, and When Not to. Bull. Math. Biol. 2019, 81, 3722–3731. [Google Scholar] [CrossRef] [Green Version]
  37. Chamseddine, I.M.; Rejniak, K.A. Hybrid modeling frameworks of tumor development and treatment. Wiley Interdiscip. Rev. Syst. Biol. Med. 2019. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Moreira, J.; Deutsch, A. Cellular Automaton Models of Tumor Development: A Critical Review. Adv. Complex Syst. 2002, 5, 247–267. [Google Scholar] [CrossRef]
  39. Lowengrub, J.S.; Frieboes, H.B.; Jin, F.; Chuang, Y.L.; Li, X.; Macklin, P.; Wise, S.M.; Cristini, V. Nonlinear modelling of cancer: Bridging the gap between cells and tumours. Nonlinearity 2010, 23, R1–R91. [Google Scholar] [CrossRef] [Green Version]
  40. Shahriyari, L. Cell dynamics in tumour environment after treatments. J. R. Soc. Interface 2017, 14, 20160977. [Google Scholar] [CrossRef]
  41. Lewin, T.D.; Maini, P.K.; Moros, E.G.; Enderling, H.; Byrne, H.M. The Evolution of Tumour Composition During Fractionated Radiotherapy: Implications for Outcome. Bull. Math. Biol. 2018, 80, 1207–1235. [Google Scholar] [CrossRef] [Green Version]
  42. Enderling, H.; Chaplain, M.A.; Anderson, A.R.; Vaidya, J.S. A mathematical model of breast cancer development, local treatment and recurrence. J. Theor. Biol. 2007, 246, 245–259. [Google Scholar] [CrossRef] [PubMed]
  43. Pérez-García, V.M.; Bogdanska, M.; Martínez-González, A.; Belmonte-Beitia, J.; Schucht, P.; Pérez-Romasanta, L.A. Delay effects in the response of low-grade gliomas to radiotherapy: A mathematical model and its therapeutical implications. Math. Med. Biol. 2015, 32, 307–329. [Google Scholar] [CrossRef] [Green Version]
  44. Kather, J.N.; Poleszczuk, J.; Suarez-Carmona, M.; Krisam, J.; Charoentong, P.; Valous, N.A.; Weis, C.A.; Tavernar, L.; Leiss, F.; Herpel, E.; et al. In Silico Modeling of Immunotherapy and Stroma-Targeting Therapies in Human Colorectal Cancer. Cancer Res. 2017, 77, 6442–6452. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Alharbi, S.A.; Rambely, A.S. A New ODE-Based Model for Tumor Cells and Immune System Competition. Mathematics 2020, 8, 1285. [Google Scholar] [CrossRef]
  46. DePillis, L.; Savage, H.; Radunskaya, A. Mathematical model of colorectal cancer with monoclonal antibody treatments. arXiv 2013, arXiv:1312.3023. [Google Scholar] [CrossRef]
  47. Sameen, S.; Barbuti, R.; Milazzo, P.; Cerone, A.; Del Re, M.; Danesi, R. Mathematical modeling of drug resistance due to KRAS mutation in colorectal cancer. J. Theor. Biol. 2016, 389, 263–273. [Google Scholar] [CrossRef]
  48. Kirshtein, A.; Akbarinejad, S.; Hao, W.; Le, T.; Su, S.; Aronow, R.A.; Shahriyari, L. Data Driven Mathematical Model of Colon Cancer Progression. J. Clin. Med. 2020, 9, 3947. [Google Scholar] [CrossRef]
  49. Beutler, B.; Greenwald, D.; Hulmes, J.D.; Chang, M.; Pan, Y.C.E.; Mathison, J.; Ulevitch, R.; Cerami, A. Identity of tumour necrosis factor and the macrophage-secreted factor cachectin. Nature 1985, 316, 552–554. [Google Scholar] [CrossRef]
  50. Grivennikov, S.I.; Greten, F.R.; Karin, M. Immunity, Inflammation, and Cancer. Cell 2010, 140, 883–899. [Google Scholar] [CrossRef] [Green Version]
  51. Mudter, J.; Neurath, M.F. IL-6 signaling in inflammatory bowel disease: Pathophysiological role and clinical relevance. Inflamm. Bowel Dis. 2007, 13, 1016–1023. [Google Scholar] [CrossRef]
  52. Waldner, M.J.; Neurath, M.F. Colitis-associated cancer: The role of T cells in tumor development. Semin. Immunopathol. 2009, 31, 249–256. [Google Scholar] [CrossRef]
  53. Terzić, J.; Grivennikov, S.; Karin, E.; Karin, M. Inflammation and Colon Cancer. Gastroenterology 2010, 138, 2101–2114. [Google Scholar] [CrossRef]
  54. Waldner, M.J.; Foersch, S.; Neurath, M.F. Interleukin-6—A key regulator of colorectal cancer development. Int. J. Biol. Sci. 2012, 8, 1248–1253. [Google Scholar] [CrossRef]
  55. Hart, A.L.; Al-Hassi, H.O.; Rigby, R.J.; Bell, S.J.; Emmanuel, A.V.; Knight, S.C.; Kamm, M.A.; Stagg, A.J. Characteristics of intestinal dendritic cells in inflammatory bowel diseases. Gastroenterology 2005, 129, 50–65. [Google Scholar] [CrossRef]
  56. Pasare, C.; Medzhitov, R. Toll pathway-dependent blockade of CD4+CD25+ T cell-mediated suppression by dendritic cells. Science 2003, 299, 1033–1036. [Google Scholar] [CrossRef] [PubMed]
  57. Aras, S.; Zaidi, M.R. TAMeless traitors: Macrophages in cancer progression and metastasis. Br. J. Cancer 2017, 117, 1583–1591. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Fan, X.; Zhang, H.; Cheng, Y.; Jiang, X.; Zhu, J.; Jin, T. Double roles of macrophages in human neuroimmune diseases and their animal models. Mediat. Inflamm. 2016, 2016, 8489251. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Deng, G. Tumor-infiltrating regulatory T cells: Origins and features. Am. J. Clin. Exp. Immunol. 2018, 7, 81–87. [Google Scholar]
  60. Iwasaki, A.; Kelsall, B.L. Freshly isolated Peyer’s patch, but not spleen, dendritic cells produce interleukin 10 and induce the differentiation of T helper type 2 cells. J. Exp. Med. 1999, 190, 229–240. [Google Scholar] [CrossRef] [Green Version]
  61. Wang, K.; Vella, A.T. Regulatory T Cells and Cancer: A Two-Sided Story. Immunol. Investig. 2016, 45, 797–812. [Google Scholar] [CrossRef]
  62. Leman, J.K.H.; Sandford, S.K.; Rhodes, J.L.; Kemp, R.A. Multiparametric analysis of colorectal cancer immune responses. World J. Gastroenterol. 2018, 24, 2995–3005. [Google Scholar] [CrossRef]
  63. Cheng, K.J.; Alshawsh, M.A.; Mejia Mohamed, E.H.; Thavagnanam, S.; Sinniah, A.; Ibrahim, Z.A. HMGB1: An overview of its versatile roles in the pathogenesis of colorectal cancer. Cell. Oncol. 2020, 43, 177–193. [Google Scholar] [CrossRef] [PubMed]
  64. Süren, D.; Yıldırım, M.; Demirpençe, Ö.; Kaya, V.; Alikanoğlu, A.S.; Bülbüller, N.; Yıldız, M.; Sezer, C. The role of high mobility group box 1 (HMGB1) in colorectal cancer. Med. Sci. Monit. Int. Med J. Exp. Clin. Res. 2014, 20, 530–537. [Google Scholar] [CrossRef] [Green Version]
  65. Guo, Z.S.; Liu, Z.; Bartlett, D.L.; Tang, D.; Lotze, M.T. Life after death: Targeting high mobility group box 1 in emergent cancer therapies. Am. J. Cancer Res. 2013, 3, 1–20. [Google Scholar]
  66. Wang, H.; Bloom, O.; Zhang, M.; Vishnubhakat, J.M.; Ombrellino, M.; Che, J.; Frazier, A.; Yang, H.; Ivanova, S.; Borovikova, L.; et al. HMG-1 as a late mediator of endotoxin lethality in mice. Science 1999, 285, 248–251. [Google Scholar] [CrossRef]
  67. Ong, S.M.; Tan, Y.C.; Beretta, O.; Jiang, D.; Yeap, W.H.; Tai, J.J.Y.; Wong, W.C.; Yang, H.; Schwarz, H.; Lim, K.H.; et al. Macrophages in human colorectal cancer are pro-inflammatory and prime T cells towards an anti-tumour type-1 inflammatory response. Eur. J. Immunol. 2012, 42, 89–100. [Google Scholar] [CrossRef] [PubMed]
  68. Darwich, L.; Coma, G.; Peña, R.; Bellido, R.; Blanco, E.J.J.; Este, J.A.; Borras, F.E.; Clotet, B.; Ruiz, L.; Rosell, A.; et al. Secretion of interferon-γ by human macrophages demonstrated at the single-cell level after costimulation with interleukin (IL)-12 plus IL-18. Immunology 2009, 126, 386–393. [Google Scholar] [CrossRef] [PubMed]
  69. Robinson, C.M.; O’Dee, D.; Hamilton, T.; Nau, G.J. Cytokines involved in interferon-gamma production by human macrophages. J. Innate Immun. 2010, 2, 56–65. [Google Scholar] [CrossRef]
  70. Zaidi, M.R.; Davis, S.; Noonan, F.P.; Graff-Cherry, C.; Hawley, T.S.; Walker, R.L.; Feigenbaum, L.; Fuchs, E.; Lyakh, L.; Young, H.A.; et al. Interferon-γ links ultraviolet radiation to melanomagenesis in mice. Nature 2011, 469, 548–553. [Google Scholar] [CrossRef] [Green Version]
  71. Liu, Y.; Zhang, P.; Li, J.; Kulkarni, A.B.; Perruche, S.; Chen, W. A critical function for TGF-beta signaling in the development of natural CD4+CD25+Foxp3+ regulatory T cells. Nat. Immunol. 2008, 9, 632–640. [Google Scholar] [CrossRef]
  72. Boyman, O.; Sprent, J. The role of interleukin-2 during homeostasis and activation of the immune system. Nat. Rev. Immunol. 2012, 12, 180–190. [Google Scholar] [CrossRef]
  73. West, N.R.; McCuaig, S.; Franchini, F.; Powrie, F. Emerging cytokine networks in colorectal cancer. Nat. Rev. Immunol. 2015, 15, 615–629. [Google Scholar] [CrossRef]
  74. Fontenot, J.D.; Rasmussen, J.P.; Gavin, M.A.; Rudensky, A.Y. A function for interleukin 2 in Foxp3-expressing regulatory T cells. Nat. Immunol. 2005, 6, 1142–1151. [Google Scholar] [CrossRef] [PubMed]
  75. Vang, K.B.; Yang, J.; Mahmud, S.A.; Burchill, M.A.; Vegoe, A.L.; Farrar, M.A. IL-2, -7, and -15, but not thymic stromal lymphopoeitin, redundantly govern CD4+Foxp3+ regulatory T cell development. J. Immunol. 2008, 181, 3285–3290. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  76. Kusume, A.; Sasahira, T.; Luo, Y.; Isobe, M.; Nakagawa, N.; Tatsumoto, N.; Fujii, K.; Ohmori, H.; Kuniyasu, H. Suppression of dendritic cells by HMGB1 is associated with lymph node metastasis of human colon cancer. Pathobiology 2009, 76, 155–162. [Google Scholar] [CrossRef]
  77. Erdman, S.E.; Poutahidis, T. Roles for Inflammation and Regulatory T Cells in Colon Cancer. Toxicol. Pathol. 2010, 38, 76–87. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  78. Ma, Y.; Shurin, G.V.; Peiyuan, Z.; Shurin, M.R. Dendritic cells in the cancer microenvironment. J. Cancer 2013, 4, 36–44. [Google Scholar] [CrossRef] [Green Version]
  79. Menetrier-Caux, C.; Montmain, G.; Dieu, M.C.; Bain, C.; Favrot, M.C.; Caux, C.; Blay, J.Y. Inhibition of the differentiation of dendritic cells from CD34+ progenitors by tumor cells: Role of interleukin-6 and macrophage colony-stimulating factor. Blood J. Am. Soc. Hematol. 1998, 92, 4778–4791. [Google Scholar]
  80. Esche, C.; Lokshin, A.; Shurin, G.V.; Gastman, B.R.; Rabinowich, H.; Watkins, S.C.; Lotze, M.T.; Shurin, M.R. Tumor’s other immune targets: Dendritic cells. J. Leukoc. Biol. 1999, 66, 336–344. [Google Scholar] [CrossRef]
  81. Sica, A.; Saccani, A.; Mantovani, A. Tumor-associated macrophages: A molecular perspective. Int. Immunopharmacol. 2002, 2, 1045–1054. [Google Scholar] [CrossRef]
  82. Sakai, Y.; Honda, M.; Fujinaga, H.; Tatsumi, I.; Mizukoshi, E.; Nakamoto, Y.; Kaneko, S. Common transcriptional signature of tumor-infiltrating mononuclear inflammatory cells and peripheral blood mononuclear cells in hepatocellular carcinoma patients. Cancer Res. 2008, 68, 10267–10279. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  83. Folkman, J.; Hochberg, M. Self-regulation of growth in three dimensions. J. Exp. Med. 1973, 138, 745–753. [Google Scholar] [CrossRef] [PubMed]
  84. Enderling, H.; Sunassee, E.; Caudell, J.J. Predicting patient-specific radiotherapy responses in head and neck cancer to personalize radiation dose fractionation. bioRxiv 2019, 630806. [Google Scholar] [CrossRef]
  85. Yu, H.; Kortylewski, M.; Pardoll, D. Crosstalk between cancer and immune cells: Role of STAT3 in the tumour microenvironment. Nat. Rev. Immunol. 2007, 7, 41–51. [Google Scholar] [CrossRef] [PubMed]
  86. Badache, A.; Hynes, N.E. Interleukin 6 inhibits proliferation and, in cooperation with an epidermal growth factor receptor autocrine loop, increases migration of T47D breast cancer cells. Cancer Res. 2001, 61, 383–391. [Google Scholar] [PubMed]
  87. Lin, M.T.; Juan, C.Y.; Chang, K.J.; Chen, W.J.; Kuo, M.L. IL-6 inhibits apoptosis and retains oxidative DNA lesions in human gastric cancer AGS cells through up-regulation of anti-apoptotic gene mcl-1. Carcinogenesis 2001, 22, 1947–1953. [Google Scholar] [CrossRef] [Green Version]
  88. Moses, H.L.; Yang, E.Y.; Pietenpol, J.A. TGF-beta stimulation and inhibition of cell proliferation: New mechanistic insights. Cell 1990, 63, 245–247. [Google Scholar] [CrossRef]
  89. Markowitz, S.D.; Roberts, A.B. Tumor suppressor activity of the TGF-beta pathway in human cancers. Cytokine Growth Factor Rev. 1996, 7, 93–102. [Google Scholar] [CrossRef]
  90. Wang, C.Y.; Eshleman, J.R.; Willson, J.K.; Markowitz, S. Both transforming growth factor-beta and substrate release are inducers of apoptosis in a human colon adenoma cell line. Cancer Res. 1995, 55, 5101–5105. [Google Scholar]
  91. Engel, M.A.; Neurath, M.F. Anticancer properties of the IL-12 family-focus on colorectal cancer. Curr. Med. Chem. 2010, 17, 3303–3308. [Google Scholar] [CrossRef]
  92. Voit, E.O.; Martens, H.A.; Omholt, S.W. 150 years of the mass action law. PLoS Comput. Biol. 2015, 11, e1004012. [Google Scholar] [CrossRef] [Green Version]
  93. Wang, Y.; Liu, C.; Liu, P.; Eisenberg, B. Field theory of reaction-diffusion: Mass action with an energetic variational approach. arXiv 2020, arXiv:2001.10149. [Google Scholar]
  94. Zoetemelk, M.; Ramzy, G.M.; Rausch, M.; Nowak-Sliwinska, P. Drug-Drug Interactions of Irinotecan, 5-Fluorouracil, Folinic Acid and Oxaliplatin and Its Activity in Colorectal Carcinoma Treatment. Molecules 2020, 25, 2614. [Google Scholar] [CrossRef] [PubMed]
  95. Wei, C.; Yang, C.; Wang, S.; Shi, D.; Zhang, C.; Lin, X.; Xiong, B. M2 macrophages confer resistance to 5-fluorouracil in colorectal cancer through the activation of CCL22/PI3K/AKT signaling. OncoTargets Ther. 2019, 12, 3051–3063. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  96. Zhang, X.; Chen, Y.; Hao, L.; Hou, A.; Chen, X.; Li, Y.; Wang, R.; Luo, P.; Ruan, Z.; Ou, J.; et al. Macrophages induce resistance to 5-fluorouracil chemotherapy in colorectal cancer through the release of putrescine. Cancer Lett. 2016, 381, 305–313. [Google Scholar] [CrossRef] [PubMed]
  97. Malesci, A.; Bianchi, P.; Celesti, G.; Basso, G.; Marchesi, F.; Grizzi, F.; Di Caro, G.; Cavalleri, T.; Rimassa, L.; Palmqvist, R.; et al. Tumor-associated macrophages and response to 5-fluorouracil adjuvant therapy in stage III colorectal cancer. OncoImmunology 2017, 6, e1342918. [Google Scholar] [CrossRef]
  98. Machover, D.; Goldschmidt, E.; Chollet, P.; Metzger, G.; Zittoun, J.; Marquet, J.; Vandenbulcke, J.M.; Misset, J.L.; Schwarzenberg, L.; Fourtillan, J.B. Treatment of advanced colorectal and gastric adenocarcinomas with 5-fluorouracil and high-dose folinic acid. J. Clin. Oncol. 1986, 4, 685–696. [Google Scholar] [CrossRef]
  99. Mini, E.; Trave, F.; Rustum, Y.M.; Bertino, J.R. Enhancement of the antitumor effects of 5-fluorouracil by folinic acid. Pharmacol. Ther. 1990, 47, 1–19. [Google Scholar] [CrossRef]
  100. Danenberg, P.V.; Gustavsson, B.; Johnston, P.; Lindberg, P.; Moser, R.; Odin, E.; Peters, G.J.; Petrelli, N. Folates as adjuvants to anticancer agents: Chemical rationale and mechanism of action. Crit. Rev. Oncol. 2016, 106, 118–131. [Google Scholar] [CrossRef] [Green Version]
  101. Kciuk, M.; Marciniak, B.; Kontek, R. Irinotecan—Still an Important Player in Cancer Chemotherapy: A Comprehensive Overview. Int. J. Mol. Sci. 2020, 21, 4919. [Google Scholar] [CrossRef]
  102. Rosner, G.; Panetta, J.; Innocenti, F.; Ratain, M. Pharmacogenetic Pathway Analysis of Irinotecan. Clin. Pharmacol. Ther. 2008, 84, 393–402. [Google Scholar] [CrossRef]
  103. de Man, F.M.; Goey, A.K.L.; van Schaik, R.H.N.; Mathijssen, R.H.J.; Bins, S. Individualization of Irinotecan Treatment: A Review of Pharmacokinetics, Pharmacodynamics, and Pharmacogenetics. Clin. Pharmacokinet. 2018, 57, 1229–1254. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  104. De Almeida, C.V.; Zamame, J.A.; Romagnoli, G.G.; Rodrigues, C.P.; Magalhães, M.B.; Amedei, A.; Kaneno, R. Treatment of colon cancer cells with 5-fluorouracil can improve the effectiveness of RNA-transfected antitumor dendritic cell vaccine. Oncol. Rep. 2017, 38, 561–568. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  105. Bailly, C.; Thuru, X.; Quesnel, B. Combined cytotoxic chemotherapy and immunotherapy of cancer: Modern times. NAR Cancer 2020, 2. [Google Scholar] [CrossRef] [Green Version]
  106. Maeda, K.; Hazama, S.; Tokuno, K.; Kan, S.; Maeda, Y.; Watanabe, Y.; Kamei, R.; Shindo, Y.; Maeda, N.; Yoshimura, K.; et al. Impact of chemotherapy for colorectal cancer on regulatory T-cells and tumor immunity. Anticancer. Res. 2011, 31, 4569–4574. [Google Scholar]
  107. Focaccetti, C.; Bruno, A.; Magnani, E.; Bartolini, D.; Principi, E.; Dallaglio, K.; Bucci, E.O.; Finzi, G.; Sessa, F.; Noonan, D.M.; et al. Effects of 5-Fluorouracil on Morphology, Cell Cycle, Proliferation, Apoptosis, Autophagy and ROS Production in Endothelial Cells and Cardiomyocytes. PLoS ONE 2015, 10, e0115686. [Google Scholar] [CrossRef]
  108. Zhang, N.; Yin, Y.; Xu, S.J.; Chen, W.S. 5-Fluorouracil: Mechanisms of Resistance and Reversal Strategies. Molecules 2008, 13, 1551–1569. [Google Scholar] [CrossRef] [Green Version]
  109. Newman, A.; Steen, C.; Liu, C.; Gentles, A.; Chaudhuri, A.; Scherer, F.; Khodadoust, M.; Esfahani, M.; Luca, B.; Steiner, D.; et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat. Biotechnol. 2019, 37, 773–782. [Google Scholar] [CrossRef]
  110. Le, T.; Aronow, R.A.; Kirshtein, A.; Shahriyari, L. A review of digital cytometry methods: Estimating the relative abundance of cell types in a bulk of cells. Briefings Bioinform. 2020. [Google Scholar] [CrossRef]
  111. Su, S.; Akbarinejad, S.; Shahriyari, L. Immune classification of clear cell renal cell carcinoma. Sci. Rep. 2021, 11, 4338. [Google Scholar] [CrossRef]
  112. Le, T.; Su, S.; Shahriyari, L. Immune classification of osteosarcoma. Math. Biosci. Eng. 2021, 18, 1879–1897. [Google Scholar] [CrossRef] [PubMed]
  113. Goldman, M.J.; Craft, B.; Hastie, M.; Repečka, K.; McDade, F.; Kamath, A.; Banerjee, A.; Luo, Y.; Rogers, D.; Brooks, A.N.; et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat. Biotechnol. 2020, 38, 675–678. [Google Scholar] [CrossRef]
  114. Kim, Y.; Friedman, A. Interaction of Tumor with Its Micro-environment: A Mathematical Model. Bull. Math. Biol. 2010, 72, 1029–1068. [Google Scholar] [CrossRef] [PubMed]
  115. Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nat. Methods 2020, 17, 261–272. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  116. FOLFIRI Chemotherapy. Available online: https://www.cancerresearchuk.org/about-cancer/cancer-in-general/treatment/cancer-drugs/drugs/folfiri (accessed on 27 March 2021).
  117. Karagiannis, G.; Hao, W.; Lin, G. Calibrations and validations of biological models with an application on the renal fibrosis. Int. J. Numer. Methods Biomed. Eng. 2020, 36, e3329. [Google Scholar] [CrossRef]
  118. Seefeld, S.; Stockwell, W.R. First-order sensitivity analysis of models with time-dependent parameters: An application to PAN and ozone. Atmos. Environ. 1999, 33, 2941–2953. [Google Scholar] [CrossRef]
  119. Yang, I.H. Uncertainty and sensitivity analysis of time-dependent effects in concrete structures. Eng. Struct. 2007, 29, 1366–1374. [Google Scholar] [CrossRef]
  120. Zi, Z. Sensitivity analysis approaches applied to systems biology models. IET Syst. Biol. 2011, 5, 336–346. [Google Scholar] [CrossRef]
  121. Dagogo-Jack, I.; Shaw, A.T. Tumour heterogeneity and resistance to cancer therapies. Nat. Rev. Clin. Oncol. 2018, 15, 81–94. [Google Scholar] [CrossRef] [PubMed]
  122. Zhu, X.; Tian, X.; Ji, L.; Zhang, X.; Cao, Y.; Shen, C.; Hu, Y.; Wong, J.W.H.; Fang, J.Y.; Hong, J.; et al. A tumor microenvironment-specific gene expression signature predicts chemotherapy resistance in colorectal cancer patients. NPJ Precis. Oncol. 2021, 5, 7. [Google Scholar] [CrossRef] [PubMed]
  123. Werner, H.M.J.; Mills, G.B.; Ram, P.T. Cancer Systems Biology: A peek into the future of patient care? Nat. Rev. Clin. Oncol. 2014, 11, 167–176. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  124. Bekisz, S.; Geris, L. Cancer modeling: From mechanistic to data-driven approaches, and from fundamental insights to clinical applications. J. Comput. Sci. 2020, 46, 101198. [Google Scholar] [CrossRef]
  125. Galon, J.; Costes, A.; Sanchez-Cabo, F.; Kirilovsky, A.; Mlecnik, B.; Lagorce-Pagès, C.; Tosolini, M.; Camus, M.; Berger, A.; Wind, P.; et al. Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science 2006, 313, 1960–1964. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  126. Taube, J.M.; Klein, A.; Brahmer, J.R.; Xu, H.; Pan, X.; Kim, J.H.; Chen, L.; Pardoll, D.M.; Topalian, S.L.; Anders, R.A. Association of PD-1, PD-1 Ligands, and Other Features of the Tumor Immune Microenvironment with Response to Anti–PD-1 Therapy. Clin. Cancer Res. 2014, 20, 5064–5074. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  127. Tredan, O.; Galmarini, C.M.; Patel, K.; Tannock, I.F. Drug Resistance and the Solid Tumor Microenvironment. JNCI J. Natl. Cancer Inst. 2007, 99, 1441–1454. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  128. Saltz, L.B.; Cox, J.V.; Blanke, C.; Rosen, L.S.; Fehrenbacher, L.; Moore, M.J.; Maroun, J.A.; Ackland, S.P.; Locker, P.K.; Pirotta, N.; et al. Irinotecan plus Fluorouracil and Leucovorin for Metastatic Colorectal Cancer. N. Engl. J. Med. 2000, 343, 905–914. [Google Scholar] [CrossRef] [PubMed]
  129. Xu, T.; Lu, J.; An, H. The relative change in regulatory T cells/T helper lymphocytes ratio as parameter for prediction of therapy efficacy in metastatic colorectal cancer patients. Oncotarget 2017, 8, 109079–109093. [Google Scholar] [CrossRef] [Green Version]
  130. Roselli, M.; Formica, V.; Cereda, V.; Jochems, C.; Richards, J.; Grenga, I.; Orlandi, A.; Ferroni, P.; Guadagni, F.; Schlom, J. The association of clinical outcome and peripheral T-cell subsets in metastatic colorectal cancer patients receiving first-line FOLFIRI plus bevacizumab therapy. OncoImmunology 2016, 5, e1188243. [Google Scholar] [CrossRef] [Green Version]
  131. Voss, A.; Voss, J. A fast numerical algorithm for the estimation of diffusion model parameters. J. Math. Psychol. 2008, 52, 1–9. [Google Scholar] [CrossRef]
  132. Parra-Rojas, C.; Hernandez-Vargas, E.A. PDEparams: Parameter fitting toolbox for partial differential equations in python. Bioinformatics 2019, 1–2. [Google Scholar] [CrossRef]
  133. Vyshemirsky, V.; Girolami, M. BioBayes: A software package for Bayesian inference in systems biology. Bioinformatics 2008, 24, 1933–1934. [Google Scholar] [CrossRef] [Green Version]
  134. Xun, X.; Cao, J.; Mallick, B.; Maity, A.; Carroll, R.J. Parameter Estimation of Partial Differential Equation Models. J. Am. Stat. Assoc. 2013, 108, 1009–1020. [Google Scholar] [CrossRef] [Green Version]
  135. Diasio, R.B.; Harris, B.E. Clinical Pharmacology of 5-Fluorouracil. Clin. Pharmacokinet. 1989, 16, 215–237. [Google Scholar] [CrossRef] [PubMed]
  136. Leucovorin Calcium Injection Label. Available online: https://www.accessdata.fda.gov/drugsatfda_docs/label/2012/040347s010lbl.pdf (accessed on 29 March 2021).
  137. Greiner, P.; Zittoun, J.; Marquet, J.; Cheron, J. Pharmacokinetics of (−)−folinic acid after oral and intravenous administration of the racemate. Br. J. Clin. Pharmacol. 1989, 28, 289–295. [Google Scholar] [CrossRef] [PubMed]
  138. Macarulla, T.; Cervantes, A.; Tabernero, J.; Roselló, S.; Van Cutsem, E.; Tejpar, S.; Prenen, H.; Martinelli, E.; Troiani, T.; Laffranchi, B.; et al. Phase I study of FOLFIRI plus pimasertib as second-line treatment for KRAS-mutated metastatic colorectal cancer. Br. J. Cancer 2015, 112, 1874–1881. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  139. Mathijssen, R.H.; van Alphen, R.J.; Verweij, J.; Loos, W.J.; Nooter, K.; Stoter, G.; Sparreboom, A. Clinical Pharmacokinetics and Metabolism of Irinotecan (CPT-11). Clin. Cancer Res. 2001, 7, 2182–2194. [Google Scholar] [PubMed]
Figure 1. Interaction network of FOLFIRI drugs. Interaction network among the key players in tumor microenvironment and FOLFIRI drugs.
Figure 1. Interaction network of FOLFIRI drugs. Interaction network among the key players in tumor microenvironment and FOLFIRI drugs.
Cancers 13 02632 g001
Figure 2. Model immune cells fractions. Clusters are obtain based on variations in 22 immune cell types of colonic tumor.
Figure 2. Model immune cells fractions. Clusters are obtain based on variations in 22 immune cell types of colonic tumor.
Cancers 13 02632 g002
Figure 3. Whole system dynamics with initial conditions as median drug doses from patient data and steady-state values for the cells. The different colour lines depict the different dynamics of each cluster. ‘Total cells’ is the sum of all immune, cancer and necrotic cells excluding cytokines.
Figure 3. Whole system dynamics with initial conditions as median drug doses from patient data and steady-state values for the cells. The different colour lines depict the different dynamics of each cluster. ‘Total cells’ is the sum of all immune, cancer and necrotic cells excluding cytokines.
Cancers 13 02632 g003
Figure 4. Varying combinations of the FOLFIRI drugs. In (A), we use the median dosage for 5-FU as the initial infusion dose, keeping the values of Leucovorin and Irinotecan zero. In (B), we add the median dose of Leucovorin to that of 5-FU, keeping Irinotecan dose at zero. Similarly in (C), 5-FU and Irinotecan are infused at the median dosage, while keeping Leucovorin at zero. (D) is the model results obtained from infusing all three drugs at median dosages. Median dosages are given in Table 3.
Figure 4. Varying combinations of the FOLFIRI drugs. In (A), we use the median dosage for 5-FU as the initial infusion dose, keeping the values of Leucovorin and Irinotecan zero. In (B), we add the median dose of Leucovorin to that of 5-FU, keeping Irinotecan dose at zero. Similarly in (C), 5-FU and Irinotecan are infused at the median dosage, while keeping Leucovorin at zero. (D) is the model results obtained from infusing all three drugs at median dosages. Median dosages are given in Table 3.
Cancers 13 02632 g004
Figure 5. Varying treatment start-time. In the corresponding (AD), the treatment is given after 1, 3, 5 and 7 years respectively. Seven years is the approximate time for the cancer to reach the steady state.
Figure 5. Varying treatment start-time. In the corresponding (AD), the treatment is given after 1, 3, 5 and 7 years respectively. Seven years is the approximate time for the cancer to reach the steady state.
Cancers 13 02632 g005
Figure 6. ROC Curves of the model. (A,B) are generated using predicted cancer cells number in the first day of follow up and the last day of follow up of the patients, and AUC values are 0.53 and 0.41, respectively.
Figure 6. ROC Curves of the model. (A,B) are generated using predicted cancer cells number in the first day of follow up and the last day of follow up of the patients, and AUC values are 0.53 and 0.41, respectively.
Cancers 13 02632 g006
Figure 7. Cancer cells dynamics with varying doses for individual patients in each cluster. (A) are obtained by varying the Leucovorin dose while keeping the 5-FU dose constant. (B) are obtained by varying 5-FU and keeping the Leucovorin dose at zero. (A1,B1), (A2,B2), (A3,B3), and (A5,B5) respectively show the results of a tumor-free patient in the clusters 1, 2, 3 and 5. Note, the patient in cluster 3 (A3,B3) is tumor-free at the first follow-up day, but is with-tumor at the last follow-up day. (A4,B4) show the results of a with-tumor patient in the cluster 5.
Figure 7. Cancer cells dynamics with varying doses for individual patients in each cluster. (A) are obtained by varying the Leucovorin dose while keeping the 5-FU dose constant. (B) are obtained by varying 5-FU and keeping the Leucovorin dose at zero. (A1,B1), (A2,B2), (A3,B3), and (A5,B5) respectively show the results of a tumor-free patient in the clusters 1, 2, 3 and 5. Note, the patient in cluster 3 (A3,B3) is tumor-free at the first follow-up day, but is with-tumor at the last follow-up day. (A4,B4) show the results of a with-tumor patient in the cluster 5.
Cancers 13 02632 g007
Figure 8. Sensitivity analysis. The first and second columns of (A) respectively present the results of non-dimensional sensitivity of cancer cell density and total cell density. (B) shows the relative sensitivity of the same quantities. Each row of plots shows the most sensitive parameters for each cluster of patients.
Figure 8. Sensitivity analysis. The first and second columns of (A) respectively present the results of non-dimensional sensitivity of cancer cell density and total cell density. (B) shows the relative sensitivity of the same quantities. Each row of plots shows the most sensitive parameters for each cluster of patients.
Cancers 13 02632 g008
Figure 9. Cancer and total cells as a function of δ C 5 f u L V . (A1A5) depict number of cancer cells and (B1B5) depict number of total cells, at different time points.
Figure 9. Cancer and total cells as a function of δ C 5 f u L V . (A1A5) depict number of cancer cells and (B1B5) depict number of total cells, at different time points.
Cancers 13 02632 g009
Figure 10. Cancer and total cells as a function of δ 5 f u M . (A1A5) depict cancer cells and (B1B5) depict Total cells, at different time points.
Figure 10. Cancer and total cells as a function of δ 5 f u M . (A1A5) depict cancer cells and (B1B5) depict Total cells, at different time points.
Cancers 13 02632 g010
Table 1. Model’s Variables. Names and descriptions of variables used in the model.
Table 1. Model’s Variables. Names and descriptions of variables used in the model.
VariableNameDescription
T N Naive T-cells
T h Helper T-cells
T C Cytotoxic cellsincludes CD8+ T-cells and NK cells
T r Regulatory T-cells
D n Naive dendritic cells
DActivated dendritic cellsantigen presenting cells
MMacrophages
CCancer cells
NNectrotic cells
HHMGB1
μ 1 Carcinogenic cytokinesincludes effects of IL-6, IL-17, IL-21 and IL-22
μ 2 Immunosuppresive agentsincludes effects of IL-10 and CCL20
I γ IFN- γ
G β TGF- β
5-FUFluorouracil
IrIrinotecan
LVLeucovorin
Table 2. The smallest tumor initial conditions in each cluster. Values of initial conditions for the dimensionless system are derived from the patients with the smallest tumor size.
Table 2. The smallest tumor initial conditions in each cluster. Values of initial conditions for the dimensionless system are derived from the patients with the smallest tumor size.
Cluster T N / T N T h / T h T C / T C T r / T r D N / D N D / D M / M
13.66 × 10 2 4.90 × 10 2 9.67 × 10 2 2.70 × 10 2 6.41 × 10 2 1.23 × 10 5 2.64 × 10 2
22.63 × 10 2 2.81 × 10 2 3.25 × 10 2 1.09 × 10 2 4.37 × 10 2 5.95 × 10 2 2.76 × 10 2
31.38 × 10 1 2.20 × 10 1 8.59 × 10 2 6.26 × 10 2 1.53 × 10 1 2.61 × 10 1 1.56 × 10 1
41.58 × 10 1 1.71 × 10 2 6.35 × 10 2 9.72 × 10 2 6.20 × 10 1 4.68 × 10 1 1.01 × 10 1
51.36 × 10 2 5.24 × 10 2 3.27 × 10 2 3.78 × 10 2 4.72 × 10 5 9.11 × 10 3 1.55 × 10 2
C / C N / N μ 1 / μ 1 μ 2 / μ 2 H / H I γ / I γ G β / G β
14.64 × 10 2 1.55 × 10 2 4.97 × 10 1 5.12 × 10 1 1.473.892.55 × 10 1
22.38 × 10 2 1.46 × 10 3 7.58 × 10 1 1.79 × 10 1 6.04 × 10 1 9.38 × 10 1 5.57 × 10 1
31.83 × 10 1 07.51 × 10 1 7.14 × 10 1 1.034.90 × 10 1 3.87
41.22 × 10 1 04.14 × 10 1 5.77 × 101.46 × 1002.56 × 10
52.69 × 10 2 7.03 × 10 3 4.59 × 10 1 2.301.184.08 × 10 1 3.46 × 10 1
Table 3. Drug dosages in mg.
Table 3. Drug dosages in mg.
FluorouracilLeucovorinIrinotecan
median770725300
min59875208
number of cycles121212
cycle length141414
Table 4. Number of patients in the validation data.
Table 4. Number of patients in the validation data.
Tumor FreeWith Tumor
Early follow up day238
Last follow up day2412
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Budithi, A.; Su, S.; Kirshtein, A.; Shahriyari, L. Data Driven Mathematical Model of FOLFIRI Treatment for Colon Cancer. Cancers 2021, 13, 2632. https://doi.org/10.3390/cancers13112632

AMA Style

Budithi A, Su S, Kirshtein A, Shahriyari L. Data Driven Mathematical Model of FOLFIRI Treatment for Colon Cancer. Cancers. 2021; 13(11):2632. https://doi.org/10.3390/cancers13112632

Chicago/Turabian Style

Budithi, Aparajita, Sumeyye Su, Arkadz Kirshtein, and Leili Shahriyari. 2021. "Data Driven Mathematical Model of FOLFIRI Treatment for Colon Cancer" Cancers 13, no. 11: 2632. https://doi.org/10.3390/cancers13112632

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop