Reaction-Diffusion on a Spatial Mathematical Model of Cancer Immunotherapy with Effector Cells and IL-2 Compounds’ Interactions

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.,ese 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.


Introduction
Cancer is one of the malignant diseases triggered by gene mutations on the cells in the tissues. e 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][2][3]. e 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. ere are some enzymes that can be used as the indicators of the mutation, e.g., p53, pRb, and EBNA1; see [1][2][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]. e 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]. e 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. e 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. is 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. e evolution profiles of the infection in a tissue are interesting to study. erefore, following the results in [11], we extend the immunotherapy model in [6][7][8] by adding the reactiondiffusion 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.

The Model Formulation
We separate the cell populations into three parts, which are the cancer cells (T), the effector cells (E), and the IL-2 compounds (I L ), as functions of time τ and position x. e diffusion terms are expressed by the second-order partial derivatives with respect to x.
In our model, we use constant parameters that have the following meanings. e 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 c and μ −1 1 , respectively. e 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 p 1 , p 2 , and p 3 , respectively. Parameters g 1 , g 2 , and g 3 , 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, s 1 represents the adoptive cellular immunotherapy (ACI) to the effector cells such as LAK or TIL cells. e cancer growth is assumed to follow the logistic model where the constant birth rate is r 2 with the carrying capacity b − 1 . Parameters μ 3 and s 3 , 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 Ω × (0, ∞), where T, E, and I L are the cancer cells, effector cells, and IL-2 compounds, respectively, and the boundary conditions are (zE/zx) � (zT/zx) � (zI L /zx) � 0 on zΩ × (0, ∞). e domain Ω is bounded in R, and zΩ is the domain domain. We adopt the formulation in [12] for the initial conditions of system (1), i.e., where E o , T o , and I L o show the natural growth rate of the effector cells, cancer cells, and the IL-2 compounds and l is the thickness of the epithelial layer. Parameters d 1 , d 2 , and d 3 , which show the random motility coefficients of the effector cells, the cancer cells, and the IL-2 compounds, are assumed to be positive. e initial conditions of the effector cell, the cancer cell, and the IL-2 compound concentrations are expressed by using the notations E 0 (x), T 0 (x), 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 Ω × (0, ∞), where the variables are and the parameters are c � (cT 0 /E 0 t s ), μ 1 � (μ 1 /t s ), p 1 � (p 1 /t s ), g 1 � (g 1 /I L 0 ), s 1 � (s 1 /t s E 0 ), r 2 � (r 2 /t s ), International Journal of Differential Equations . e transformation is motivated by the results in [13]. e initial condition of system (3) is and the boundary conditions are on (0, ∞).

The Cancer-Free Equilibrium
e 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., s 1 � 0 and s 3 � 0, and the others are for the treatment cases, i.e., s 1 > 0 and s 3 � 0, s 1 � 0 and s 3 > 0, and s 1 > 0 and s 3 > 0.
e proof of the theorem is given by using the linearization of system (3) near the equilibrium point e 0 i with i � 1, 2, 3, 4. See [14] for the method.

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. e numerical data that we used in this paper are based on the clinical data from [6,12,[15][16][17]. Based on [6], the value of the parameter c is 0 ≤ c ≤ 0.2778, μ 1 � 0.167, .001, and μ 3 � 55.55556. Following [12,[15][16][17], we use Following [15], we use t s � 0.18. 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).
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.

The Dynamics of the System for Nontreatment Cases
We consider the dynamics of system (3) for s 1 � s 3 � 0. 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 eorem 1, we found that the cancer-free equilibrium e 0 1 is unstable. For the numerical simulation, based on [6], we use c � 0.25( ∼ 0.045 day − 1 ) that shows antigenicity for cancer, and the random motility d 2 is equal to 0.0000199( ∼ 3.6 × 10 − 6 cm 2 day − 1 ). Parameter p 2 � 0.5555556 shows the decay rate of the cancer cells by the effector cells, and the diffusion coefficients d 1 � 0.001 and d 3 � 0.01. We show the results in Figure 2.
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. e 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 eorem 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 d 2 as d 2 � 0.00048( ∼ 10 − 9 cm 2 s − 1 ), then the evolution profiles of each concentration are changed. In this case, the solutions of system (3) cannot be isolated into a certain profile.
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 t, and we found that the solutions converge to a stable limit cycle. e 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.

The Immunotherapy Case by Using ACI and Cytokine Therapy
In this section, we consider the case that s 1 > 0 and s 3 > 0 that represent the immunotherapy using ACI and cytokine therapy. In this simulation, we use s 1 � 0.0035 and s 3 � 0.2. e 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. ese results are due to the fact that, for the immunotherapy case, if the conditions in eorem 1 have been satisfied, we have the stable cancer-free equilibrium.   Distance into tissue (i) Figure 5: Evolution profiles of the effector cell (blue), the cancer cell (green), and IL-2 (red) concentrations for the cancer cell random motility coefficient  Figures 6 and 7. We vary the values of the diffusion parameter that shows the motility of the cancer cells, i.e., d 2 � 0.0000048, d 2 � 0.0000199, d 2 � 0.00025, and d 2 � 0.00048. 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 c � 0.25 for the higher-antigenicity case where the profiles are shown in Figure 6. For the lower-antigenicity case, we use c � 0.000475, and the profile is presented in Figure 7.
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 c � 0.000475. In this case, the cancer cells' concentration becomes higher and covers the large part of the tissue; see Figure 7. e 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 d 2 , 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.

Concluding Remarks
e 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. ere are two important parameters in our system that represent the antigenicity for cancer and the motility of the cancer cells. e 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. e limit cycle behaviour is usually caused by the Hopf bifurcation when the parameter value is varied. is 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.
e analysis of this case is one of the open problems for our system. e 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 e 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.  Figure 7: e spatial-temporal evolutions of the cancer cells' diffusion coefficients for t � 25, t � 50, and t � 100 with antigenicity parameter c � 0.000475. e top row is for d 2 � 0.0000048, the second row is for d 2 � 0.0000199, the third row is for d 2 � 0.00025, and the last row is for d 2 � 0.00048.

Conflicts of Interest
International Journal of Differential Equations 9 in UGM and the Cancer Modelling Team, UGM, for the discussions during the research.