Abstract
Immunotherapy is one of the future treatments applicable in most cases of cancer including malignant cancer. Malignant cancer usually prevents some genes, e.g., p53 and pRb, from controlling the activation of the cell division and the cell apoptosis. In this paper, we consider the interactions among the cancer cell population, the effector cell population that is a part of the immune system, and cytokines that can be used to stimulate the effector cells called the IL-2 compounds. These interactions depend on both time and spatial position of the cells in the tissue. Mathematically, the spatial movement of the cells is represented by the diffusion terms. We provide an analytical study for the constant equilibria of the reaction-diffusion system describing the above interactions, which show the initial behaviour of the tissue, and we conduct numerical simulation that shows the dynamics along the tissue that represent the immunotherapy effects. In this case, we also consider the steady-state conditions of the system that show the long-time behaviour of these interactions.
1. Introduction
Cancer is one of the malignant diseases triggered by gene mutations on the cells in the tissues. The gene mutations affect the cell cycle, DNA damages, and some anomalies on the cells. Some anomalies’ behaviour of the enzymes in the cells due to the DNA damage caused by cancer has been shown in [1–3].
The mathematical model that shows the dynamics of cancer on the tissues has been given in [4, 5]. In these papers, the authors consider the cervical cancer case and show the boundaries, in the state space and in the parameter space, which separate two types of solutions. One solution goes to an equilibrium point, and the other one grows up to infinity that represents the cancer cells’ metastases. When the precancerous growth rate parameter and the saturation parameter of the cancer cells are varied, the boundary on the state space is adjusted; see [4].
In the subcellular level, the cancer infections are mainly caused by the gene mutation represented by the shifting concentration of the enzymes. There are some enzymes that can be used as the indicators of the mutation, e.g., p53, pRb, and EBNA1; see [1–3]. For the breast cancer case, the gene mutations triggered by the enzyme reactions can cause anomalies on the cell cycle and cell repair regulations; see [2, 3].
The immunotherapy model of cancer that involves the interaction between the cancer cells, the effector cells, and the IL-2 compounds was motivated by Kirschner et al. [6, 7]. The basic model, the interactions, and the global stability of the equilibria have been introduced in [6, 7]. In [8], the authors added the periodic perturbation term in the system that shows the periodic stimulation of the effector cells by the IL-2 compounds and studied the bifurcations of the system due to the variation of the perturbation parameter.
The other immunotherapy method for cancer is done by using oncolytic virus; see [9, 10]. Virus is injected into the cancer site to reduce the DNA replication error while increasing the apoptosis rate of the cancer cells. This type of immunotherapy also plays an important role for stimulating the effector cells and cytokines to prevent the growth of cancer cells.
In the real situation, the growth direction of the cancer cells depends on the weakest parts of the tissue. The evolution profiles of the infection in a tissue are interesting to study. Therefore, following the results in [11], we extend the immunotherapy model in [6–8] by adding the reaction-diffusion terms on each component. Unlike the model in [11], we use different coordinate nondimensional transformations to prevent the Michaelis–Menten constants on the system to understand the role of these parameters in the dynamics of the system.
Our model is a three-dimensional system of partial differential equations, and the analysis is focused on the study of the steady state and the cycle of infection of the disease that is represented by the limit cycle of the system. We apply the Runge–Kutta method to determine the solutions of the system numerically. It is important to understand the behaviour of the system.
2. The Model Formulation
We separate the cell populations into three parts, which are the cancer cells (), the effector cells (), and the IL-2 compounds (), as functions of time and position . The diffusion terms are expressed by the second-order partial derivatives with respect to .
In our model, we use constant parameters that have the following meanings. The antigenicity for the cancer cells that measures the ability of the immune system to recognize cancer via the non-self-protein antigens and the average of the natural lifespan of the effector cells are represented by parameters and , respectively. The maximum proliferation rate of the effector cells by the IL-2 compounds, the maximum degradation rates of the cancer cells that have an interaction with the effector cells, and the maximum production rates of the IL-2 compounds by the effector cells are represented by parameters , , and , respectively. Parameters , , and , respectively, denote the Michaelis–Menten constants that represent the proliferation kinetics of the effector cells, the degradation kinetics of the cancer cells, and the growth of the IL-2 compounds caused by the interaction between the cancer cells and the effector cells. Moreover, represents the adoptive cellular immunotherapy (ACI) to the effector cells such as LAK or TIL cells.
The cancer growth is assumed to follow the logistic model where the constant birth rate is with the carrying capacity . Parameters and , respectively, show the decay rate of the IL-2 compounds and the cytokine therapy to increase IL-2 compounds.
Our model is formulated as a three-dimensional system of PDE as follows:in , where , and are the cancer cells, effector cells, and IL-2 compounds, respectively, and the boundary conditions are on . The domain is bounded in , and is the domain domain. We adopt the formulation in [12] for the initial conditions of system (1), i.e.,where , and show the natural growth rate of the effector cells, cancer cells, and the IL-2 compounds and is the thickness of the epithelial layer. Parameters , , and , which show the random motility coefficients of the effector cells, the cancer cells, and the IL-2 compounds, are assumed to be positive. The initial conditions of the effector cell, the cancer cell, and the IL-2 compound concentrations are expressed by using the notations , , and .
We apply the nondimensional transformation for the variables and parameters of system (1) to study the interactions between these variables. We have the transformed system as the following:in , where the variables areand the parameters are , , , , , , , , , , , , , , and , . The transformation is motivated by the results in [13].
The initial condition of system (3) isand the boundary conditions areon .
3. The Cancer-Free Equilibrium
The analysis of the cancer-free equilibrium is important to determine the conditions of the patients to be cured from cancer. For the numerical simulation, we focus on the dynamics of the solutions near the cancer-free equilibrium.
System (3) has a cancer-free equilibrium solution, i.e.,whose stability depends on the treatment parameters which are applied to the system. We have four cases of the stability for the cancer-free equilibrium: one is for the nontreatment case, i.e., and , and the others are for the treatment cases, i.e., and , and , and and .
Theorem 1. (1)If and , then the equilibrium is unstable.(2)If and , then the equilibrium is locally asymptotically stable.(3)If and , then the equilibrium is unstable.(4)If and , then the equilibriumis locally asymptotically stable.
The proof of the theorem is given by using the linearization of system (3) near the equilibrium point with . See [14] for the method.
4. The Initial Position of the Effector Cells, Cancer Cells, and IL-2 Compounds on the Tissue
In the beginning of cancer invasions, the cancer cells enter the tissue via the epithelial layer. In this case, the effector cells and the IL-2 compounds, which are parts of the immune system, will respond to attack the cancer cells before they enter the body via the outer epithelial layer.
The numerical data that we used in this paper are based on the clinical data from [6, 12, 15–17]. Based on [6], the value of the parameter is , , , , , , , , and . Following [12, 15–17], we use , , , , and . Following [15], we use . It implies that the time unit on system (1) is equal to 5.555556 times of the time unit of system (3).
In Figure 1, we show the initial profile of the effector cells, cancer cells, and IL-2 compounds. In Figure 1, we found that the cancer cells are found on the inner epithelial layer, which are between 0 and 0.2 cm; see Figure 1(b). At the same time, the effector cells and the IL-2 compounds, as parts of the immune system, are found on the outer epithelial layer which are between 0.2 and 1 cm; see Figures 1(a) and 1(c).
(a)
(b)
(c)
In the following sections, we perform some numerical simulations to study the dynamics of system (1) that represents spread of the cancer cells in a tissue. We employ the Runge–Kutta method and adopt the method in [18] to perform numerical simulations. In this case, we will simulate the change of the positions of the cancer cells, the effector cells, and the IL-2 compounds several times and show the steady-state patterns of the solutions.
5. The Dynamics of the System for Nontreatment Cases
We consider the dynamics of system (3) for . In this case, we will show the natural patterns of the effector cells and IL-2 compounds in response to the appearance of the cancer cells.
According to Theorem 1, we found that the cancer-free equilibrium is unstable. For the numerical simulation, based on [6], we use that shows antigenicity for cancer, and the random motility is equal to . Parameter shows the decay rate of the cancer cells by the effector cells, and the diffusion coefficients and . We show the results in Figure 2.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
In Figure 2, we show that the effector cell and the IL-2 compound concentrations give a response for the appearance of the cancer cells, which is to isolate the activities of the cancer cells. The increase of the effector cells’ concentration will stimulate the increase of the IL-2 compounds’ concentration. By the fact that the cancer-free equilibrium is unstable (see Theorem 1), it is predicted that the cancer cells’ concentration cannot be vanished from the system. However, it is isolated to a certain profile on the tissue.
In Figure 3, we show that if we adjust the random motility as , then the evolution profiles of each concentration are changed. In this case, the solutions of system (3) cannot be isolated into a certain profile.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
By the results in Figures 2 and 3, we found that the random motility of the cancer cells can change the evolution profile of each concentration on the tissue. In Figure 4, we show the solution of system (3) with respect to time , and we found that the solutions converge to a stable limit cycle. The stable limit cycle in this case represents that the concentrations of the effector cells, cancer cells, and the IL-2 compounds always fluctuate, but they are bounded to a certain periodic cycle.
(a)
(b)
6. The Immunotherapy Case by Using ACI and Cytokine Therapy
In this section, we consider the case that and that represent the immunotherapy using ACI and cytokine therapy. In this simulation, we use and .
The evolution profile of the effector cell, cancer cell, and the IL-2 compound concentrations is shown in Figure 5. By using the immunotherapy, we found that the cancer cells’ concentration is decreasing and then vanishes after sometime. These results are due to the fact that, for the immunotherapy case, if the conditions in Theorem 1 have been satisfied, we have the stable cancer-free equilibrium.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
7. The Spatial-Temporal Evolutions of the Cancer Cells Based on the Variation of the Diffusion Parameter
The evolution profile of the cancer cells’ concentrations for the nontreatment case is shown in Figures 6 and 7. We vary the values of the diffusion parameter that shows the motility of the cancer cells, i.e., , , , and . In this case, we will compare the evolution profile of the higher and the lower antigenicity for the cancer cells. We use the antigenicity parameter for the higher-antigenicity case where the profiles are shown in Figure 6. For the lower-antigenicity case, we use , and the profile is presented in Figure 7.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
From the results in Figure 6, we found that the cancer cells cannot vanish after sometime. However, in the long period of time, the cancer cells’ concentration in the tissue is isolated into an intermediate level.
Different situations occur in system (3) when we have the lower antigenicity for the cancer cells, that is, for . In this case, the cancer cells’ concentration becomes higher and covers the large part of the tissue; see Figure 7. The situations are due to the fact that most of the receptors of the immune system, which is represented by the effector cells and the IL-2 compounds, cannot detect activities of cancer.
From Figures 6 and 7, we found the fact that if the motility of the cancer cells, which is represented by the diffusion parameter , increases, then the spread of the cancer cells in the tissue becomes faster. In the bottom right of Figure 7, we show that most of the tissue is covered by the cancer cells.
8. Concluding Remarks
The spread of the cancer cells on the tissue does not only depend on the time but also on the position of the cancer cells in the tissue. There are two important parameters in our system that represent the antigenicity for cancer and the motility of the cancer cells. The higher motility of the cancer cells implies that the spread of the cancer cells in the tissue becomes faster. In this case, the lower antigenicity causes higher concentration of the cancer cells to spread faster in the tissue. For the higher antigenicity, although the cancer cells cover most of the tissue, the cells’ concentration is on the intermediate level.
One of the important phenomena in our system is the appearance of the stable limit cycle when we choose a certain value of the diffusion parameter of the cancer cells. The limit cycle behaviour is usually caused by the Hopf bifurcation when the parameter value is varied. This bifurcation is one of the entry points to the possibility for the system to have chaotic behaviour that represents the metastases of the cancer cells. The analysis of this case is one of the open problems for our system.
The other open problem is that, in the real situation, the cancer cells can grow and spread to other tissues and organs via blood vessels. In this case, we can approach the model using the moving boundary condition that determines the spread of the cancer cells in the body.
Data Availability
The authors do not use the data from any repositories to provide the results of the research. For the numerical simulations, the value of parameters are adopted from the reference [6] and by the assumptions.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was partially supported by the Deputy for Research and Strengthening Development of the Ministry of Research and Technology, Indonesia (the National Research and Innovation Agency), which has provided research funding through Penelitian Dasar Unggulan Perguruan Tinggi (PDUPT) with contract agreement numbers 6/AMD/E1/KP.PTNBH/2020 and 2773/UN1.DITLIT/DIT-LIT/PT/2020. The authors would also like to thank some colleagues in UGM and the Cancer Modelling Team, UGM, for the discussions during the research.