1 Introduction

Colon cancer typically originates as polyps in the intestinal wall and affects the large intestine, Accounting for approximately 9.4 percent of all cancer-related deaths in 2020. Symptoms of the disease, most commonly digestive problems, may not appear until the advanced stages [1, 10, 11, 19]. Within the tumor’s inflammatory microenvironment, the cellular component includes host immune cells such as neutrophils, dendritic cells, mast cells, monocytes/macrophages, natural killer cells, and types of T cells. It has been shown that CD8\(^+\) T cells and effector CD4\(^+\) T cells may have antitumor properties [3]. In [2], it is demonstrated that tumor rejection relies on the CD8\(^+\) T cell subset. The results indicated that CD4\(^+\) CD25\(^+\) immunoregulatory cells played a role in the progression of leukemias, myeloma, and sarcomas [4]. In instances where tumor-infiltrating lymphocytes (TILs) have enhanced patient outcomes, T lymphocytes are acknowledged as the primary agents of anti-tumor immune responses. A specific group of CD4\(^+\) T cells, known as CD4\(^+\) CD25\(^+\) regulatory T cells (Treg), can gather in the tumor microenvironment and inhibit tumor-specific T-cell reactions, consequently impeding tumor elimination [5]. The cellular makeup of the tumor microenvironment and the presence of CD8\(^+\), CD4\(^+\), and Regulatory T (T\(_{reg}\)) cells are discussed in [6]. Effector T cells need a secondary activation step from tissue dendritic cells in order to attain full proficiency for effector functions, proliferation, and survival. The stimulation of effector T cells by dendritic cells within tissues has evolved as a mechanism to guarantee that T cells are activated to their maximum capacity only at the location of the active infection [7]. The existence of tumor emboli in lymphovascular and perineural structures is regarded as an initial stage of metastatic invasion. There is a correlation between the presence of an immune response within the tumor and the lack of tumor emboli. Specifically, colorectal cancer without tumor emboli was linked to a heightened infiltration of immune cells and an upregulation of mRNA expression of adaptive T-helper 1 (TH1) effector T-cell markers [8].

Dendritic cells (DCs) are crucial components of the immune system, responsible for presenting antigens to T cells and initiating immune responses. In the context of colon cancer, DCs have been found to have a complex role, exhibiting both tumor-inhibiting and tumor-promoting properties. It is evident that dendritic cells play a central role in the immune system due to their capability to regulate both immune tolerance and immunity. Therefore, DCs are a crucial focus in endeavors to induce therapeutic immunity against cancer [9].

The interactions between the immune system, including DCs, and colon tumor cells are complex and not fully understood. Understanding these interactions is crucial to comprehend the mechanisms of tumor immune surveillance and evasion, as well as metastasis, and to develop new approaches for treating colon cancer. Dendritic cells have been investigated as a potential immunotherapy for colon cancer. DCs vaccines, which entail extracting DCs from a patient’s blood, loading them with tumor antigens, and re-infusing them into the patient, have exhibited promise in clinical trials [21].

The dynamics of biological phenomena, including diseases, can be better understood using time-dependent differential equations. Rules and equations governing biochemical reactions enable a detailed examination of the interactions between cancer cells and the body’s immune system components, most of which are in the form of ordinary differential equations (ODEs). The use of ODEs to describe different aspects of cancer, including tumor dynamics, the tumor microenvironment, interactions with the immune system, and predicting its progression, has been in practice since the early 1990 s [57, 60, 63, 65, 67, 69, 70]. A tumor-immune interactions model, focusing on the role of natural killer (NK) and CD8+ T cells in tumor surveillance, to understand the dynamics of immune-mediated tumor rejection is presented in [57]. Various ways of the circulatory system and rheological properties of blood have been described in a mathematical way [60]. Vitale et al. [65, 70] have described a model for B-CLL that assumes competition between cancer cells and the immune response, wherein a single combined population of T lymphocytes represents the entire immune response. Examining the nonlinear dynamics of immunogenic tumors, the model developed by Kuznetsov [63] showcases oscillatory growth patterns in tumors, along with dormancy and ’creeping through’−a phase where the tumor remains minuscule for an extended duration before rapidly expanding to a critical size. These mathematical models attribute the cyclical behavior of the tumor directly to its interaction with the immune system. Various mathematical models, including those proposed by [12,13,14, 22], have been used to investigate and predict the course of colon cancer.

A data-based mathematical model of colon cancer in the form of an ODE system has been developed in [3] to investigate the development of immune system components such as T cells, DCs, and macrophages. The authors divided the data of cancer patients into five clusters based on their immune system, generated five virtual patients using the data from each cluster, and built their tumor growth model. In this paper, we will use the data from three different clusters to construct a mathematical model of colon cancer. Fractional-order differential equation models are another category of mathematical models that offer advantages regarding relevant data, providing information about memory and hereditary properties, and being non-local in nature dependent on initial values. Due to these properties, fractional-order models are preferred over classical models in describing real phenomena [23,24,25,26,27]. The viscoelastic properties of various biological phenomena such as the liver, brain, heart, valves, special cells, and calcium dynamics can often be described using fractional-order differential equations [35,36,37,38,39,40,41,42,43, 59, 62, 64, 66]. Recently, fractional models have been effectively utilized to describe the behavior of B-CLL cancer and COVID-19, with promising results published [28,29,30]. In this paper, we introduce a fractional-order model to understand the interactions between cancer cells and the immune system, as well as to determine the effective factors in tumor progression and metastasis. Introduced by Raissi et al. [15,16,17], Physics-Informed Neural Networks (PINNs) are a type of neural network designed to solve nonlinear partial differential equations (PDEs). Training neural networks can enable the solution of supervised learning tasks while respecting the physical laws described by general PDEs. There are two types of PINNs: the forward problem, which involves a data-driven solution of differential equations, and inverse problems, which entail the data-driven discovery of dynamical systems. In the second type of PINN, neural networks are used to identify suitable parameters for the dynamical system. In the second type of PINN, neural networks are used to identify suitable parameters for the dynamical system. Disease-Informed Neural Networks (DINNs) are a state-of-the-art method inspired by PINNs that incorporates differential equations in a mathematical model of the disease [18]. They have demonstrated the effectiveness of DINNs by using the SIR-model to learn the dynamics of several families of infectious diseases and forecast their progression. Epidemiology-Informed Neural Network (EINN) was introduced in [20] to solve the asymptomatic SIR model, which shows the proportion of asymptomatic infected individuals to the total number of infected individuals. Furthermore, to solve stochastic differential equations [32, 33], differential equations of fractional order and integro-differential equations [31], the PINNs have been extended. The open-source library DeepXDE is written in Python and implements the PINN algorithm [34]. In this paper, we use a mathematical model based on ordinary differential equations of fractional order to better understand the complexities of colon cancer and the impact of the body’s immune system in dealing with it. In order to avoid excessive complexity, we only model some of the key players and interaction networks involved in colon cancer, as determined through our review of articles on colon cancer progression, despite the involvement of numerous cell types and molecules. Based on the network depicted in Fig. 1, the resulting mathematical model has been developed. After building the model, to estimate the parameters of the model and obtain the solution of the model, the PINN method will be developed by the Fractional-Cancer-Informed Neural Networks (FCINN)-based deep learning as an effective method for solving system of fractional differential equations. The organization of the paper is as follows. In Sect. 2, we develop a mathematical model of fractional order for the immune system response to colon cancer. The analysis of the model is presented in Sect. 3. We extend the PINN to FCINN for parameter estimation and solving the model in Sect. 4. The discussion and results are presented in Sect. 5. In conclusion, a summary is presented in Sect. 6.

2 The mathematical model

Regarding fractional derivatives, there are a couple of commonly used definitions: Grunwald-Letnikov, Caputo, and Riemann-Liouville, etc. [68].

Definition 2.1

[24] The Caputo fractional derivative of order \(\alpha\) of g(t), \(t>0\) is defined as

$$ D^{\alpha }g(t)=\frac{1}{\Gamma (n-\alpha )} \int _{0}^{t} (t- \tau )^{n-\alpha -1}f^{n}(\tau ) d\tau ,$$

where \(\Gamma (\cdot )\) is the Gamma function, \(\alpha \in (n-1,n)\), \(n \in \mathbb {Z}^{+}\).

The Caputo derivative Laplace transform for \(n-1\le \alpha < n\) is described as:

$$\begin{aligned} \mathcal {L} \{ D^{\alpha } g(t);s\}=s^{\alpha } G (s)-\sum _{j=1}^{n-1} s^{\alpha -i-1} g^{(j)}(0), \end{aligned}$$

where \(G(s) = \mathcal {L} \{g(t)\}\). If \(g^{(j)}(0) = 0\) for \(j = 1, 2, \dots , n-1\), then the Laplace transform of the Caputo derivative can be simplified to \(\mathcal {L} \{ D^{\alpha } g(t);s\}=s^\alpha G (s)\).

Fig. 1
figure 1

Cellular network

Full size image

The fractional-order model considers the population of colon cancer cells as well as the populations of three compartments of the immune system: dendritic cells, cytotoxic T cells, and T-helper cells. The variable time (t) is assumed to represent the duration of time in days. To build the model, we consider T(t), D(t), L(t), and H(t) as the populations of total colon cancer cells, total dendritic cells, cytotoxic T cells (CD\(8^+\) T-cells), and total helper T cells (CD\(4^+\) T-cells), respectively. The populations of these cells are typically measured as concentrations in units of cells per cubic centimeter (cm\(^3\)). Figure 1 provides a graphic representation of the cellular network dynamics that describe the interactions between the colon tumor and the three populations of the immune system. This figure is an important visual aid that helps to illustrate the behavior of the model and the mechanisms underlying the immune response to colon cancer. The proposed cellular network and mathematical model in this paper are based on a set of assumptions (A1–A6) that are valid medical facts. These assumptions provide a foundation for the model and help to ensure that it accurately reflects the behavior of the system being studied.

  1. A1)

    It is well established that the growth and proliferation of cells follow a logistic pattern [45, 46].

  2. A2)

    Demonstrated in some studies that dendritic cells have the ability to kill tumor cells [9, 48].

  3. A3)

    Studies have demonstrated that cytotoxic T cells have the ability to inhibit tumor cells [5, 8].

Based on the assumptions (A1–A3) previously described, the population dynamics of colon tumor cells can be modeled as follows:

$$\begin{aligned} \frac{d^\alpha T(t)}{dt^\alpha } =pT(1-qT)-d_{LT}LT-d_{DT}DT, \end{aligned}$$
(1)

where \(\alpha \in (0,1)\) represents the fractional order, p represents the tumor cell proliferation rate, \(q^{-1}\) is the total capacity of tumor cells, and \(d_{LT}\) and \(d_{DT}\) represent the inhibiting rates of tumor cells by CD\(8^+\) T-cells and DCs, respectively. These inhibitions are modeled according to the mass action law.

  1. A4)

    Demonstrated in some studies that cancer cells can activate dendritic cells, and DCs have the ability to lyse tumor cells [47, 48].

The dynamics of the DCs population can be described by the following equation:

$$\begin{aligned} \frac{d^\alpha D(t)}{dt^\alpha } =s_{1}+a_{TD}\mathcal{F}(T)D-d_{TD}\mathcal{G}(T)D-d_DD, \end{aligned}$$
(2)

where \(s_D\) is the production rate of DCs, \(a_{TD}\) and \(d_{TD}\) are the activation and inactivation rates of dendritic cells by tumor cells, respectively, and \(d_D\) represents the death rate of dendritic cells. In Eq. 2, \(\mathcal{F}(T)=\frac{T^2}{m_1+T^2}\) is a Hill function and \(\mathcal{G}(T)=\frac{T}{m_2+T}\) is a Michaelis-Menten function.

  1. A5)

    The recruitment of CD\(8^+\) T-cells is affected by newly activated CD\(4^+\) T-cells [51, 52].

  2. A6)

    CD\(8^+\) T-cells activation increases with DCs [48, 49].

  3. A7)

    Activated CD\(8^+\) T-cells have been shown to enhance the inactivation of colon cancer cells by expressing high levels of cytokines such as IFN-gamma and FasL, which enable them to kill the cancer cells [48, 50, 51, 53].

Based on assumptions (A5–A7), we can derive the following equation to describe the dynamics of cytotoxic T cells:

$$\begin{aligned} \frac{d^\alpha L(t)}{dt^\alpha } =a_{HL}\mathcal{M}(T)H+k_{TD}\mathcal{N}(T)D-d_{TL}TL-d_LL, \end{aligned}$$
(3)

where \(\mathcal{M}(H)=\frac{H}{m_3+H}\) and \(\mathcal{N}(T)=\frac{T}{m_4+T}\) terms are in the form of Michaelis-Menten functions. The parameter \(a_{HL}\) represents the activation rate of CD\(8^+\) T-cells by newly activated CD\(4^+\) T-cells, \(d_{TL}\) is the clearing rate of cancer cells by CD\(8^+\) T-cells, and \(d_L\) represents the death rate of CD\(8^+\) T-cells.

  1. A8)

    Dendritic cells can activate helper T cells through antigen presentation [51].

The dynamics of CD\(4^+\) T-cells can be described by the following equation:

$$\begin{aligned} \frac{d^\alpha H(t)}{dt^\alpha } =s_{H}+a_{DH}\mathcal{W}(D)L-d_HH. \end{aligned}$$
(4)

The recruitment of CD\(4^+\) T-cells by dendritic cells is modeled according to the Michaelis-Menten function with rate \(a_{DH}\) and \(\mathcal{W}(D)=\frac{D}{m_5+D}\). Helper T cells can be produced at a constant rate of \(s_H\) and are lost at a rate of \(d_H\). The resulting dynamics of the interactions between colon cancer cells and the immune system (equations 1-4) can be described by the following system of fractional differential equations:

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{d^\alpha T(t)}{dt^\alpha } =pT(1-qT)-d_{LT}LT-d_{DT}DT,\\ \frac{d^\alpha D(t)}{dt^\alpha } =s_{1}+a_{TD}\mathcal{F}(T)D-d_{TD}\mathcal{G}(T)D-d_DD,\\ \frac{d^\alpha L(t)}{dt^\alpha } =a_{HT}\mathcal{M}(T)H+k_{TD}\mathcal{N}(T)D-d_{TL}TL-d_LL,\\ \frac{d^\alpha H(t)}{dt^\alpha } =s_{H}+a_{DH}\mathcal{W}(D)L-d_HH, \end{array}\right. } \end{aligned}$$
(5)

and the initial conditions are \(T(0)\ge 0,~D(0)\ge 0,~L(0)\ge 0,~H(0)\ge 0\).

Here, \(\alpha \in (0,1)\) represents the fractional order. In system (5), the first equation represents the dynamics of colon tumor cells, the second equation describes the DCs dynamics, the third equation describes the dynamics of cytotoxic T-cells, and the last equation describes the dynamics of helper T-cells. The parameters and functions used in the equations are defined as previously described.

3 Analysis of the model

3.1 The non-negativity of the solutions in the model

Ensuring the non-negativity of solutions in biological models is necessary.

Lemma 3.1

(Generalized Mean Value Theorem [54]) Assuming that g(t) is a continuous function on [a, b] and the fractional derivative \(\frac{d^\alpha g(t)}{dt^\alpha }\) is continuous for real values of a, b, and \(\alpha \in (0, 1]\), then we have:

$$\begin{aligned} \frac{g(t)-g(a)}{(t-a)^\alpha }\Gamma ({\alpha })=\frac{d^\alpha }{dt^\alpha }g(\eta ), ~~0\le \eta \le t,~~ \forall t \in (a, b]. \end{aligned}$$

Lemma 3.2

If all initial values of system (5) are non-negative, then their solutions for all positive times t are also non-negative.

Proof

Let \(\mathcal{D}=\lbrace (T, D, L, H)\in \mathbb {R}^4| T\ge 0, D\ge 0, L\ge 0, H\ge 0\rbrace\) be a positively invariant region. Then, on the hyperplanes within the region \(\mathcal{D}\), we have

$$\begin{aligned} \frac{d^\alpha }{dt^\alpha }T(t)|_{T=0}=0,&~~~~~~~~ \frac{d^\alpha }{dt^\alpha }D(t) |_{D=0} = s_D\ge 0,\\ \frac{d^\alpha }{dt^\alpha }L(t) |_{L=0} =0,&~~~~~~~~ \frac{d^\alpha }{dt^\alpha }H(t) |_{H=0} = s_{H}+a_{DH}\mathcal{W}(D)L\ge 0.\\ \end{aligned}$$

If \(\lbrace T(0), D(0), L(0), H(0) \rbrace \in \mathcal{D}\) then according to Lemma 3.1 we have

$$\begin{aligned} \frac{d^\alpha }{dt^\alpha }T(t) \ge 0,~\frac{d^\alpha }{dt^\alpha }D(t) \ge 0,~\frac{d^\alpha }{dt^\alpha }L(t) \ge 0, \frac{d^\alpha }{dt^\alpha }H(t) \ge 0, \end{aligned}$$

therefore, the solutions T(t), D(t), L(t), H(t) are increasing functions and they will remain on the hyperplanes \(\mathcal{D}\) forever. Thus, if all initial values of system (5) are non-negative, then its solutions for all positive times t are also non-negative. \(\square\)

3.2 Stability of the equilibrium points

It is necessary to study the long-term solutions of the model to gain insight into the dynamics of colon cancer. To achieve this goal, we calculate the equilibrium points of the model with and without cancer, which are the values for which the system changes are zero [29]. Cancer-free equilibrium of the model (5) is given by \(E_0=(0, \frac{s_D}{d_D}, 0, \frac{s_H}{d_H})\), other equilibrium points \((T^*, D^*, L^*, H^*)\) can be obtained numerically by solving following algebraic system:

$$\begin{aligned} \frac{d^\alpha }{dt^\alpha }T(t)=\frac{d^\alpha }{dt^\alpha }D(t)=\frac{d^\alpha }{dt^\alpha }L(t)=\frac{d^\alpha }{dt^\alpha }H(t)=0. \end{aligned}$$

Due to biological considerations, we only consider non-negative equilibrium points. For any equilibrium point, we linearize the system to analyze its behavior.

$$\begin{aligned} \frac{d^\alpha }{dt^\alpha }T(t)=\;&\Big (p-2pqT^*-d_{LT}L^*-d_{DT}D^*\Big )T-d_{LT}T^*L-d_{DT}T^*D, \nonumber \\ \frac{d^\alpha }{dt^\alpha }D(t)=\;&\Big (a_{TD}D^*\frac{\partial \mathcal{F}}{\partial T}(T^*)-d_{TD}D^*\frac{\partial \mathcal{G}}{\partial T}(T^*)\Big )T+\Big (a_{TD}\mathcal{F}(T^*)-d_{TD}\mathcal{G}(T^*)-d_D\Big )D, \nonumber \\ \frac{d^\alpha }{dt^\alpha }L(t)=\;&\Big (a_{TH}H^*\frac{\partial \mathcal{M} }{\partial T}(T^*)+k_{TD}D^*\frac{\partial \mathcal{N}}{\partial T}(T^*)-d_{TL}L^*\Big )T \nonumber \\&+a_{TH}\mathcal{M}(T^*)H+k_{TD}\mathcal{N}(T^*)D-\Big (d_{TL}T^*+d_L\Big )L, \\ \frac{d^\alpha }{dt^\alpha }H(t)=&a_{HD}H^*D\frac{\partial \mathcal{W}}{\partial D}(D^*)+\Big (a_{HD}\mathcal{W}(H^*)-d_H\Big )H.\nonumber \end{aligned}$$
(6)

By applying the Laplace transform on the sides of (6), we have

$$\begin{aligned} s^\alpha \mathcal{T}(s)-s^{\alpha -1}T(0)=\;&\Big (p-2pqT^*-d_{LT}L^*-d_{DT}D^*\Big )\mathcal{T}-d_{LT}T^*\mathcal{L}-d_{DT}T^*\mathcal{D}, \nonumber \\ s^\alpha \mathcal{D}(s)-s^{\alpha -1}D(0)=\;&\Big (a_{TD}D^*\frac{\partial \mathcal{F}}{\partial T}(T^*)-d_{TD}D^*\frac{\partial \mathcal{G}}{\partial T}(T^*)\Big )\mathcal{T}+\Big (a_{TD}\mathcal{F}(T^*)-d_{TD}\mathcal{G}(T^*)-d_D\Big )\mathcal{D}, \nonumber \\ s^\alpha \mathcal{L}(s)-s^{\alpha -1}L(0)=\;&\Big (a_{TH}H^*\frac{\partial \mathcal{M} }{\partial T}(T^*)+k_{TD}D^*\frac{\partial \mathcal{N}}{\partial T}(T^*)-d_{TL}L^*\Big )\mathcal{T} \nonumber \\&+a_{TH}\mathcal{M}(T^*)\mathcal{H}+k_{TD}\mathcal{N}(T^*)\mathcal{D}-\Big (d_{TL}T^*+d_L\Big )\mathcal{L}, \nonumber \\ s^\alpha \mathcal{H}(s)-s^{\alpha -1}H(0)=&a_{HD}H^*\frac{\partial \mathcal{W}}{\partial D}(D^*)\mathcal{D}+\Big (a_{HD}\mathcal{W}(H^*)-d_H\Big )\mathcal{H}, \end{aligned}$$
(7)

where \(\mathcal{T}(s), \mathcal{D}(s), \mathcal{L}(s), \mathcal{H}(s)\) are Laplace transforms of T(t), D(t), L(t), H(t), respectively. System (7) can be written in the following matrix form:

$$\begin{aligned} \Xi (s)* \begin{pmatrix} \mathcal{T}(s)\\ \mathcal{D}(s) \\ \mathcal{L}(s)\\ \mathcal{H}(s) \end{pmatrix}= \begin{pmatrix} z_1(s)\\ z_2(s)\\ z_3(s)\\ z_4(s) \end{pmatrix}, \end{aligned}$$
(8)

where \(z_1(s)=s^{\alpha -1}T(0), z_2(s)=s^{\alpha -1}D(0), z_3(s)=s^{\alpha -1}L(0), z_4(s)=s^{\alpha -1}H(0)\) and

$$\begin{aligned}\Xi (s)= \begin{pmatrix} s^{\alpha }+\xi _1 &{} \xi _2 &{} \xi _3 &{} 0\\ \xi _4 &{} s^{\alpha }+\xi _5 &{} 0&{} 0 \\ \xi _6 &{}\xi _7 &{} s^{\alpha }+\xi _8 &{} \xi _9\\ 0 &{} \xi _{10} &{} 0 &{} s^{\alpha }+\xi _{11} \end{pmatrix}. \end{aligned}$$

\(\Xi (s)\) is the characteristic matrix for system (5) at equilibrium point \((T^*, D^*, L^*, H^*)\) and

$$\begin{aligned} \xi _1 = p-2pqT^*-d_{LT}L^*-d_{DT}D^*,~~ \xi _2 =d_{LT}T^*,~~ \xi _3 = d_{DT}T^*,\\ \xi _4 = a_{TD}D^*\frac{\partial \mathcal{F}}{\partial T}(T^*)-d_{TD}D^*\frac{\partial \mathcal{G}}{\partial T}(T^*),~~ \xi _5 = a_{TD}\mathcal{F}(T^*)-d_{TD}\mathcal{G}(T^*)-d_D,\\ \xi _6 =a_{TH}H^*\frac{\partial \mathcal{M} }{\partial T}(T^*)+k_{TD}D^*\frac{\partial \mathcal{N}}{\partial T}(T^*)-d_{TL}L^*,~~\xi _7 = k_{TD}\mathcal{N}(T^*),\\ \xi _8 = -d_{TL}T^*-d_L,~~\xi _9 =a_{TH}\mathcal{M}(T^*),~~\xi _{10} = a_{HD}H^*\frac{\partial \mathcal{W}}{\partial D}(D^*),~~\xi _{11}=a_{HD}\mathcal{W}(H^*)-d_H. \end{aligned}$$

Eigenvalues of the characteristic matrix \(\Xi (s)\) at cancer-free equilibrium point \(E_0\) are

$$\begin{aligned} \lambda _1=p-\frac{d_{DT}s_D}{d_D},~~~\lambda _2=-d_D,~~~\lambda _3=-d_L,~~~\lambda _4=a_{HD}\mathcal{W}(\frac{s_H}{d_H})-d_H. \end{aligned}$$

The equilibrium point \(E_0\) is stable if \(p<\frac{d_{DT}s_D}{d_D}\) and \(a_{HD}\mathcal{W}(\frac{s_H}{d_H})<d_H.\)

Lemma 3.3

The equilibrium point \(E=(T^*, D^*, L^*, H^*)\) is locally asymptotically stable if \(\beta _0>0, \beta _1>0, \beta _3>0, \beta _1\beta _2\beta _3>\beta _1^2+\beta _3^2\beta _0.\)

Proof

The characteristic equation at \(E=(T^*, D^*, L^*, H^*)\) is defined by

$$\begin{aligned} s^{4\alpha }+\beta _3 s^{3\alpha }+\beta _2 s^{2\alpha }+\beta _1 s^{\alpha }+\beta _0 =0, \end{aligned}$$
(9)

\(\beta _3=\xi _1+\xi _2+\xi _3+\xi _4,~\beta _2=\xi _1\xi _2+\xi _3\xi _4+\xi _3+\xi _4,~\beta _1=\xi _3\xi _4(\xi _1+\xi _2)+\xi _1\xi _2(\xi _3+\xi _4)\), \(\beta _0=\xi _1\xi _2\xi _3\xi _4.\) Based on Routh-Hurwitz criterion, E is an equilibrium point of the system that is stable in a local neighborhood if \(\beta _0>0, \beta _1>0, \beta _3>0, \beta _1\beta _2\beta _3>\beta _1^2+\beta _3^2\beta _0.\) \(\square\)

4 Fractional-cancer-informed neural networks

We have developed the FCINN, a deep learning framework that utilizes the fractional differential system (5) to estimate model parameters and compute numerical solutions. The system captures the dynamics of the interactions between colon cancer and the immune system. A schematic representation of the cancer-informed neural networks is shown in Fig. 2.

Fig. 2
figure 2

Schematic representation of FCINN

Full size image

The FCINN formula uses a combined deep neural network-dynamic system technique and a loss function to simulate the mathematical model. The loss function in this method consists of two parts. The loss function consists of two parts: The first one measures the discrepancy between the neural network output and the available data, while the second one accounts for the residuals of the fractional-ODE. The network parameters are optimized for learning data by minimizing the sum of these two parts. At the same time, by minimizing the residual error of the dynamic system, we can also deduce the parameters of the mathematical model. Since automatic differentiation technology cannot be used for fractional operators directly, the FCINN structure for models with fractional derivatives is typically different from classical models. Therefore, we first discretize the fractional derivative using finite difference before performing FCINN. This process allows us to apply automatic differentiation to the resulting discretized system and obtain the necessary gradients for optimizing the network parameters. To simplify and condense the notation, we write the system of equations (5) as follows:

$$\begin{aligned} \frac{d^{\alpha }}{dt^{\alpha }}{\mathbf{x (t)}}=\textbf{f}(\mathbf{{x}}, t,\varvec{\mu }),~~~\mathbf{{x}}(t_0)=\mathbf{{x}}_0, \end{aligned}$$
(10)

where \(\textbf{x}=[T, D, L, H]^T\) represents the concentration of dependent variables, and the model parameters that need to be determined are

$$\begin{aligned} \varvec{\mu }=(p, q, d_{LT}, d_{DT}, s_D, a_{TD}, d_{TD}, d_D, a_{HT}, k_{TD}, d_{TL}, d_L, s_H, a_{DH}, d_H). \end{aligned}$$

In the deep neural network part of FCINN, we use a feed-forward neural network that can be represented as a function of the input vector \(\textbf{t}\), K hidden layers, and output:

$$\begin{aligned} NN(t; \theta )= \sigma ({W}_K\sigma (...\sigma (W_2\sigma (W_1t+b_1)+b_2)...)+b_K, \end{aligned}$$

where \(\varvec{\theta }\) includes weight matrices \(\textbf{W}_i\) and bias vectors \(\textbf{b}_i\) that are the trainable parameters for the i-th layer, and \(\sigma\) is a nonlinear activation function. We consider the sigmoid activation function \(\sigma (t) = \frac{1}{1+e^t}\). The mean squared error loss function for the neural network is defined as

$$\begin{aligned} MSE_{NN}=\frac{1}{M}\sum _{i=1}^{5}\sum _{j=1}^M||\mathbf{{NN}}_i(t_j; \varvec{\theta })-\textbf{x}_{ij}||_2^2, \end{aligned}$$

where \(\textbf{NN}_i(t_j;\varvec{\theta })=[\tilde{T}(t_j), \tilde{D}(t_j), \tilde{L}(t_j), \tilde{H}(t_j)]^T\) represents the target variables estimated by the deep neural network at time points \(t_j\), \(j=1,\dots ,M\). The training dataset denoted by \(\textbf{x}_{ij}\) is obtained from the given dataset in [3]. Based on FCINN, the outputs of the neural network (\(\tilde{T}, \tilde{D}, \tilde{L}, \tilde{H}\)) must satisfy the residual of the cancer-informed dynamical system. Therefore, we consider the following equation

$${\mathbf{Eq}}({\mathbf{t}},{\varvec{\mu }}): = \widehat{{\mathbf{x}}}({\text{t}}) - {\mathbf{f}}({\mathbf{x}},{\mathbf{t}}, {\varvec{\mu}} ),$$
(11)

where \({\hat{{\mathbf{x}}}}(t)\) is the discretized form of \(\frac{d^{\alpha }}{dt^{\alpha }}\textbf{x}(t)\) that is used to apply automatic differentiation. The mean squared error loss function for the cancer-informed system is given by:

$$\begin{aligned} MSE_{Eq}=\frac{1}{M}\sum _{i=1}^{5}\sum _{j=1}^M|| Eq_i(t_j, {\boldsymbol{\theta}} , {\boldsymbol{\mu }}) ||_2^2. \end{aligned}$$
(12)

The total mean squared error loss of FCINN is defined as a function of both sets of parameters \(\varvec{\theta }\) and \(\varvec{\mu }\):

$$\begin{aligned} MSE=MSE_{NN}+MSE_{Eq}. \end{aligned}$$
(13)

Through the utilization of gradient-based optimizers, we can use the loss function (13) to infer both the neural network parameters \(\varvec{\theta }\) and the unknown parameters of the model \(\varvec{\mu }\). Additionally, the numerical solutions of the model (5) can be accurately predicted by the trained neural network. The use of gradient-based optimization allows us to efficiently search the parameter space and find the optimal values that minimize the loss function and improve the accuracy of the predictions.

5 Discussion and results

5.1 Patients’ data

To estimate the model parameters and validate our approach, we utilized patient data obtained from the deconvolution of gene expression within colon tumors by Arkadz Kirshtein et al. [3]. To estimate the relative frequency of immune cells within each tumor, the authors utilized the CIBERSORTx B-mode on gene expression profiles. They then classified patients into different groups based on their immune patterns. Using the machine learning algorithm of K-means clustering, based on the analysis of the immune patterns, the authors could distinguish five distinct groups of colon cancers. For our study, we used data from three of these clusters. These data were valuable in developing and validating our approach for modeling the interactions between colon cancer and the immune system.

5.2 Parameters estimation

To estimate the model parameters, we used data related to the colon cancer cell population of three out of five patients in [3] for training and carried out the FCINN algorithm. For the parameter identification part of the PINN, we employed three hidden linear layers, each consisting of 5 neurons. On the other hand, for the compartment size prediction part, we utilized six hidden layers, each comprising 350 neurons. The activation function ”sigmoid” is utilized in each hidden layer, while the ”softplus” activation function \(\phi (x)=ln(1+e^x)\) is employed in the output layer. We chose the Adam algorithm as the optimizer for the task. The results of fitting one sample data are displayed in Fig. 3 (left), and some estimated parameters are presented in Table 1.

Fig. 3
figure 3

Data fitting for the first cluster (left), Simulations of colon cancer cells in 3 clusters (right)

Full size image
Table 1 Parameter values in the model
Full size table

As shown in Fig. 3, we obtain reasonable fits to colon cancer data using our model equations, which suggests that our approach is effective in modeling the interactions between colon cancer and the immune system. The time scale was normalized by multiplying it by 0.04 for consistency across all patients’ data.

5.3 Sensitivity analysis

Several model parameters have a significant impact on determining the model output, and even minor changes in these parameters can cause drastic changes in the predictions. To identify these parameters and their level of influence on controlling the progression of colon cancer, we will use sensitivity analysis techniques. By doing so, the level of influence of each component of the immune system included in the model on controlling the progression of colon cancer can be determined. To quantify the sensitivity and relationships between the model parameters, we utilized the Partially Ranked Correlation Coefficients (PRCCs) method. PRCC was calculated for 1000 values of each parameter generated using the Latin Hypercube Sampling (LHS) method, as described in [56]. LHS is a type of Monte Carlo sampling that allows for the simultaneous adjustment of all input parameters. This sampling method is effective when the outcome is a monotonic function of each input parameter. In our study, we only included the parameters a, \(d_{LT}\), \(s_D\), \(a_{TD}\), \(d_{TD}\), \(k_{TD}\), \(a_{DH}\), \(a_{HT}\), \(d_{DT}\), and the initial value of tumor density \(T_0\) in the sensitivity analysis. These parameters were found to be monotonically associated with the model outcomes, and their PRCC values provide valuable insights into the sensitivity and relationships between the model parameters. The sensitivity analysis results can help us identify the most critical parameters for predicting colon cancer progression.

Fig. 4
figure 4

Parameter correlation values

Full size image

5.4 Numerical simulations

Fractional-order differential equation models have advantages in predicting changes in colon tumor cells based on the characteristics of the immune system. One of these advantages is demonstrated by the ability to consider different fractional-order levels and their effects on the patterns of changes in colon tumor cells, the immune system, and their interactions. In our study, we obtained three distinct clusters of the disease by applying the FCINN method for three different values: \(\alpha =0.98\), \(\alpha =0.88\), and \(\alpha =0.78\), as described in [3]. The simulation results are presented in Fig. 3 (right).

Figure 4 shows that the parameter determining the proliferation rate of cancer cells, p, and the initial amount of cancer cells, \(T_0\), have a positive correlation with the population of colon cells, indicating that an increase in these parameters leads to an increase in the number of cancer cells. In contrast, the parameters determining the rates of proliferation (\(s_D\)) and recruitment of dendritic cells \((a_{TD})\), as well as the inhibition of cancer cells by dendritic cells (\(d_{DT}\)), have the strongest negative correlation with the population of colon cancer cells. This indicates that an increase in these parameters leads to a decrease in the number of cancer cells. Targeting the proliferation rate of cancer cells or increasing the recruitment and activity of dendritic cells may be effective strategies for treating colon cancer.

To investigate the effects of the identified influential parameters more accurately, we performed simulations of the model solutions using the FCINN method for various values of these parameters. Specifically, we obtained solutions of the model for three values of dendritic cell activation rate by tumor, \(a_{TD}=0.01, 0.03, 0.06\), which are shown in Fig. 5. The different plots in Fig. 5 show the changes in the population of cancer cells and immune cells in three simulated clusters. In all three clusters, the cellular populations reach the same steady-state values but with different patterns. For the value of \(a_{TD}=0.01\), the equilibrium value of the model is equal to

$$\begin{aligned} \textbf{e}_1=(4.5\times 10^6, 7.13\times 10^3, 3.30\times 10^3, 5.05\times 10^3). \end{aligned}$$

The order of the fractional derivative in the model affects the rate at which the solutions of the model reach the equilibrium value of the system, with higher fractional orders resulting in faster convergence. In the first cluster (\(\alpha =0.98\)), the population of cancer cells reaches its maximum value after the 330th day of tumor diagnosis, while in the second (\(\alpha =0.88\)) and third (\(\alpha =0.78\)) clusters, they reach their maximum values after the 500th and 900th days, respectively. These results demonstrate the importance of considering the fractional order of the derivative in modeling the dynamics of colon cancer and its interaction with the immune system. We conducted simulations of the behaviors of colon tumor cells and the immune system for \(a_{TD}=0.03\), and the results are shown in Fig. 6. It can be observed that the plots for the population of cancer cells show a decrease in the population of cancer cells in all three clusters as the activation rate \(a_{TD}\) increases. Solutions of the model approach the equilibrium point

$$\begin{aligned} \textbf{e}_2=(123.4, 7.62\times 10^5, 8.29\times 10^3, 6.73\times 10^3). \end{aligned}$$

However, the population of cancer cells in different clusters reaches the equilibrium point at different times. In the first cluster, the population of cancer cells decreases after 600 days from the first day of tumor diagnosis. In the second and third clusters, the population of cancer cells starts to decrease on the 1400th and 2600th days, respectively. These results provide important insights into the effects of the activation rate \(a_{TD}\) on the dynamics of colon cancer and its interaction with the immune system.

The plots in Fig. 7 show that increasing the \(a_{TD}\) parameter to 0.06 leads to faster inhibition of tumor cells, and the model solutions tend to the equilibrium point sooner. We also conducted a stability analysis of the equilibrium point \(\textbf{e}_2\) by examining the real parts of the eigenvalues of the characteristic matrix \(\Xi (s)\) at the equilibrium point. We found that all the real parts of the eigenvalues are negative, indicating that the equilibrium point is stable. This conclusion is also supported by simulating the model responses for different values of \(a_{TD}\).

Interestingly, we found that the value of \(a_{TD}=0.0237\) represents a bifurcation point for the activation rate of dendritic cells against the tumor. Specifically, for values lower than 0.0237, the population of cancer cells reaches its maximum value, indicating that the immune system is unable to control tumor growth. However, for values greater than 0.0237, tumor growth is inhibited, indicating that the immune system is effectively controlling tumor growth. This finding highlights the critical role of dendritic cells in the immune response to colon cancer and provides important insights into the dynamics of colon cancer and its interaction with the immune system. It can guide future research and treatment strategies aimed at enhancing the activation and recruitment of dendritic cells for the treatment of colon cancer.

We examined the effect of parameter \(d_{DT}\), the tumor inhibition rate by dendritic cells, on the behavior of the model solutions, and the results are presented in the plot of Fig. 8. As shown in the plot, even a slight change in this parameter has a significant impact on reducing the time required to inhibit the population of colon cells. Specifically, we considered three values of \(d_{DT}=10^{-6},10^{-7}, 10^{-8}\) in cluster three. For \(d_{DT}=10^{-8}\), the tumor starts to decrease after day 800, while this decrease occurs after day 720 and 440 for \(d_{DT}=10^{-7}\) and \(d_{DT}=10^{-6}\), respectively. In Fig. 9, it is shown that a critical production rate of the DCs population exists. Numerical solutions of the model with \(s_D\) value below this critical value, the tumor cell population grows to an equilibrium point. Conversely, solutions with \(s_D\) value above the critical value grow to a zero equilibrium.

These results highlight the critical role of DCs in the immune response to colon cancer and demonstrate the importance of considering the tumor inhibition rate by dendritic cells in modeling the dynamics of colon cancer and its interaction with the immune system. They can guide future research and treatment strategies aimed at enhancing the activity and efficacy of dendritic cells for the treatment of colon cancer. One approach involves taking dendritic cells from a patient’s blood and modifying them in the laboratory to make them more effective at recognizing and attacking cancer cells. These modified dendritic cells are then injected back into the patient’s body, where they can activate an immune response against the cancer. Another approach is to use vaccines that contain dendritic cells and tumor antigens to stimulate the immune system to attack colon cancer cells. This can be done either by injecting the vaccine directly into the tumor or by injecting it into the patient’s bloodstream. Compared to study [3], with a much smaller model and fewer equations, valuable and interesting results were obtained in this work. The effect of dendritic cells in the process of colon disease was determined more precisely.

Fig. 5
figure 5

Simulation results for \(a_{TD}=0.01\) in three clusters

Full size image
Fig. 6
figure 6

Simulation results for \(a_{TD}=0.03\) in three clusters

Full size image
Fig. 7
figure 7

Simulation results for \(a_{TD}=0.06\) in three clusters

Full size image
Fig. 8
figure 8

Simulation results for different values of \(d_{DT}=0.01\) in third cluster

Full size image
Fig. 9
figure 9

Simulation results for different values of \(s_{D}=0.01\) in third cluster

Full size image

6 Conclusion

In this study, we investigated the behavior of colon cancer cell population in long-term interactions with three cell populations: dendritic cells, cytotoxic T cells, and helper T cells, using fractional differential equations. This modeling approach has the advantage of incorporating memory and non-locality, allowing for clearer visualization of the system’s state relative to its temporal changes compared to ordinary differential equation models. Moreover, we demonstrated that different fractional orders can produce different behavioral patterns of the immune system and cancer progression, enabling us to identify three patient clusters based on their interactive behavior of colon cancer cells with the immune system.

To estimate the model parameters and solve the equations, we developed an algorithm based on the deep learning of PINN for the FCINN method, which is suitable for simulating fractional-order models. By implementing this algorithm, we obtained high-accuracy parameter estimates and detailed simulations of the model behavior over the long term. Our sensitivity analysis identified dendritic cell-related parameters as having the highest correlation with the population of cancer cells, which can guide future research and treatment strategies for colon cancer. The sensitivity analysis showed that some parameters of the model have a significant effect on the process of colon cancer (Fig. 4). The proliferation rate of cancer cells, p, and the initial amount of cancer cells, \(T_0\), have a positive correlation with the population of colon cells. The recruitment of dendritic cells (\(a_{TD}\)), as well as the inhibition of cancer cells by dendritic cells (\(d_{DT}\)), has the strongest negative correlation with the population of colon cancer cells.

Furthermore, we performed simulations using the FCINN method to confirm the effectiveness of dendritic cells in the control and treatment of colon cancer. We also determined the equilibrium points of the model and their stability and identified the bifurcation point, which corresponds to the critical value of the activation of dendritic cells by tumor cells parameter, to determine the qualitative behavior of the model. This critical amount can be used as the optimal value of the fraction of dendritic cells in the treatment of colon cancer.  The influence of the parameter, \(a_{TD}\), dendritic cell activation rate by tumor, on the changes in the population of cancer cells and the immune system in three clusters is demonstrated by the solutions of the model for three values of \(a_{TD} = 0.01, 0.03, 0.06\), as shown in Figs. 5, 6 and 7 . It can be observed that the plots for the population of cancer cells show a decrease in the population of cancer cells in all three clusters as the activation rate \(a_{TD}\) increases.  By performing bifurcation analysis, the bifurcation point for \(a_{TD}\) parameter is 0.0237. For values lower than the bifurcation point, the tumor cell population reaches its maximum value (Fig. 5).  Conversely, for values greater than the bifurcation point, the growth of tumor cells is controlled, reaching its lowest value, which represents the steady state of the system (Figs. 6 and 7). This finding highlights the critical role of dendritic cells in the immune response to colon cancer. The effect of parameter \(d_{DT}\), the tumor inhibition rate by dendritic cells, on the behavior of the model solutions is presented in the plot of Fig. 8. For \(d_{DT} = 10^8\), the tumor starts to decrease after day 800, while this decrease occurs after day 720 and 440 for \(d_{DT} = 10^7\) and \(d_{DT} = 10^6\), respectively. These results highlight the critical role of DCs in the immune response to colon cancer and demonstrate the importance of considering the tumor inhibition rate by dendritic cells in modeling the dynamics of colon cancer and its interaction with the immune system. They can guide future research and treatment strategies aimed at enhancing the activity and efficacy of dendritic cells for the treatment of colon cancer. Notably, dendritic cell-based immunotherapy has been shown to improve the overall survival of patients with metastatic colorectal cancer [55]. Therefore, our study suggests that dendritic cells may have the potential to inhibit the metastatic properties of colon cancer cells and improve patient outcomes. Overall, while dendritic cell therapy for colon cancer is still in its early stages, it shows promise as a potential new treatment option for this disease.

In summary, our study highlights the potential of using fractional differential equations to optimize the application of immunotherapy for colon cancer treatment. By incorporating memory and non-locality, and by considering different fractional orders, we were able to identify patient clusters and critical parameters for effective treatment strategies. Our findings provide valuable insights into the complex dynamics of colon cancer and its interaction with the immune system and can guide future research and treatment approaches for this disease.