# Mathematical and Numerical Modelling of Interference of Immune Cells in the Tumour Environment

Authors: Sweta Sinha; Paramjeet Singh; Mehmet Emir Koksal
Journal: Discrete Dynamics in Nature and Society
 (2023)

---

## Abstract

In this article, the behaviour of tumour growth and its interaction with the immune system have been studied using a mathematical model in the form of partial differential equations. However, the development of tumours and how they interact with the immune system make up an extremely complex and little-understood system. A new mathematical model has been proposed to gain insight into the role of immune response in the tumour microenvironment when no treatment is applied. The resulting model is a set of partial differential equations made up of four variables: the population density of tumour cells, two different types of immune cells (CD4+ helper T cells and CD8+ cytotoxic T cells), and nutrition content. Such kinds of systems also occur frequently in science and engineering. The interaction of tumour and immune cells is exemplified by predator-prey models in ecology, in which tumour cells act as prey and immune cells act as predators. The tumour-immune cell interaction is expressed via Holling’s Type-III and Beddington-DeAngelis functional responses. The combination of finite volume and finite element method is used to approximate the system numerically because these approximations are more suitable for time-dependent systems having diffusion. Finally, numerical simulations show that the methods perform well and depict the behaviour of the model.


---

## Body

## 1. Introduction

A tumour is a collection of cancer cells that grows due to uncontrolled cell proliferation. The immune system has the ability to control such growth of tumours. It can distinguish tumour cells from other cells and, to some extent, eliminate them. Such an immune system belongs to two main elements: innate and adaptive. The innate response provides a quick-to-mount general defence mechanism. This response includes natural killer cells, macrophages, and neutrophil leukocytes. Unlike the innate response, the adaptive response takes longer (up to 7 days) to activate a defence, but it targets abnormal cells more precisely once initiated. Adaptive element cells include dendritic cells, B lymphocytes (B cells), and T lymphocytes (T cells). The initiation of immune cells to eliminate tumour cells begins with the release of the tumour-specific antigen from dying tumour cells.Furthermore, immune cells like dendritic cells or macrophages help transport the antigen to the lymph node. Subsequently, T cells recognize the antigens present in the lymph node and travel to the tumour site via blood vessels. After reaching, it infiltrates the tissue to attack tumour cells, releasing more tumour-specific antigens and thus repeating the tumour immunity interaction process [1]. Figure 1 depicts the entire process of the immune system penetrating tumour tissue.Figure 1
The tumour-immune interaction cycle (this figure is taken from [1]).Many clinical studies focus on enhancing the tumour-immune interaction process by strengthening the immune system with vaccines or directly injecting T cells or cytokines (immunotherapy) [2]. Moreover, Owolabi et al. [3] emphasized the interaction between transformed epithelial cells (TECs), fibroblasts, myofibroblasts, transformed growth factor (TGF), and epithelial growth factor, which can aid in the administration of drugs to destroy tumour cells. Kim et al. [4] discussed efficient strategies to reduce the count of tumour cells. However, before developing any immunotherapies, it is necessary to comprehend the mechanism by which a growing tumour interacts with the immune system, which is presently unclear. Mathematical modelling has been a valuable tool in exploring the significance of tumour-immune interactions. The various tumour dynamics can be studied using mathematical modelling. For example, [5, 6] studied the mathematical model of tumour dynamics combined with the chemotherapeutic treatment and mathematically analyzed the model for the stability requirement. The mechanism underlying the tumour-immune interaction can be expressed using a mathematical model, whose details can be found in a review article [7].The research presented in [8] gives a detailed description of the multiple frameworks used to model the dynamics of the tumour-immune system. Many researchers have described nonspatial and spatial models to study the interaction. For instance, Matzavinos et al. [9] present a spatial model describing the penetration of lymphocytes (T cells and B cells) in solid tumours. The model examined some crucial parameters, including the inactivation or death of immune cells when targeting tumour cells, the binding rate of the immune cells with the tumour cells, and the chemotactic response of the immune cells in the presence of chemokines (a large group of proteins helpful in simulating migratory behaviour of immune cells), which is useful in assessing cancer’s dormancy. In addition, [10] investigated the interaction between the immune system and a tumour that was expanding as a result of a diffuse nutrition supply using the spatial model. Whereas in the article [11], the authors proposed a new model of tumour-immune cell interactions using a system of three coupled ordinary differential equations (nonspatial). The main focus was to examine the role of natural killer and CD8+ T cells in tumour detection. Moreover, [12, 13] modifies the nonspatial model by introducing a constant time delay in the growth rate of the immune system’s hunting cells. The one possible approach to capture the tumour-immune interactions comes from ecology that uses the concept of predator-prey models [14–16]. In such a model, the immune cells act as predators, while the tumour cells play the role of prey. The analogy may not be directly applicable to tumour-immune interaction. For example, (i) there is no direct conversion of tumour cells into immune cells, (ii) the immune cell survival does not rely on the number of tumour cells, and (iii) in mature tumours, the prediction of oscillatory behaviour is uncommon. Nevertheless, the modifications of an ecological class of models serve as a starting point for developing the tumour-immune interaction model that may help with targeted treatment. As a result, we propose a new mathematical model which studies the spatial dynamics of the interaction of predator (immune cells) with prey (tumour cells).The study of functional response is critical in understanding the interaction of immune cells and tumour cells in terms of the prey-predator approach. Various functional responses explain the interaction between the tumour-immune system. For instance, Agrawal et al. [17] have considered Hollings-II functional response [18] to explain the tumour-immune interaction. In the present study, Holling’s type-III [19] and Beddington-DeAngelis (BD) functional responses [20] are used to define the interaction. Holling’s type-III functional response motivated by Dawes and Souza [18] is expressed for tumour-immune interaction as(1)fT=aT2h2+T2,where h denotes the half-saturation constant, a is the immune cell’s maximal attacking rate, and T is the tumour cell density. The BD functional response, which is initially introduced by Beddington [21] and De Angelis et al. [22], on the other hand, is expressed as(2)fT,I=αTλ+T+γI,where I denotes immune cell density, and the positive constants α, γ, and λ represent the rate of maximum attacking rate, saturation constant, and immune cell interference parameter, respectively.Furthermore, in the immune system, T cells are an important part of the immune system in the elimination process of tumour cells. It develops in the bone marrow and matures in the thymus (a small organ which produces T cells). In the thymus, T cells differentiate into CD4+ helper, regulatory, or CD8+ cytotoxic T cells. Once stimulated by the appropriate antigen, CD8+ cytotoxic T cells kill cancer cells, regulatory T cells control immune reactions, and CD4+ helper T cells facilitate numerous cell types. Therefore, our model only considers two types of T cells (CD8+ cytotoxic T cells and CD4+ helper T cells). Both T cells act complementary in the process of elimination of tumour cells. The behaviour is explained clinically and incorporated into many mathematical models. For example, Dritschel et al. [23] examined the role of helper T cells and cytotoxic T cells when the immune system interacts with the tumour. Kirschner and Panetta proposed a system of three ODEs to explain the dynamics between the population of T cells, tumour cells, and cytokine IL-2 produced by the CD4+ helper T cells [24]. According to their model, if tumour cells have a low antigenicity, the immune system will have difficulty clearing them. The immune system can reduce the tumour size if tumour cell antigenicity is high. On the other hand, Robertson-Tessi et al. described an interaction between CD4+ helper T cells, other cells from the immune system, and tumour cells using twelve ODEs [25]. Their model captures the three E’s (elimination, equilibrium, and escape) of immunoediting tumour cells. Furthermore, the dynamics of tumour-immune interaction were studied by incorporating the immune activation delay (time delay) for effector cells, and CD4+ helper T cells [26, 27]. This time delay helped in recognizing the stability region for a tumour cell.In this paper, we study the prey-predator type model to understand the interaction of tumour cells and two types of immune cells (CD4+ helper T cells and CD8+ cytotoxic T cells). The model presented by El-Gohary [15] and Agarawal et al. [17] is the basis for the contribution in the present model. The modifications are listed below:(i)
The nonlinear system of the ordinary differential equation is extended to the nonlinear system of the partial differential equation. This will give insight into each cell’s spatial distribution and migration patterns during the interaction.(ii)
The nutrients are present in sufficient quantities during the early stages of tumour formation. Therefore, nutrient-dependent tumour growth is included. Moreover, the chemotaxis of tumour cells, directed migration of tumour cells towards the gradient of nutrient concentration, is incorporated.(iii)
The functional response is modified from Holling’s type-II to BD functional response for expressing the interaction of tumour cells and CD8+ cytotoxic T cells.(iv)
In this work, we have considered the activation and exhaustion rate of CD8+ cytotoxic T cells due to interaction with tumour cells which is expressed by Holling’s type-III functional response.(v)
The CD4+ helper T cells stimulate the proliferation of CD8+ cytotoxic T cells, which is expressed using Holling’s type-III functional response.(vi)
Instead of using a PDE solver, we computed the solution of the model variable using the classical numerical methods, the finite volume method (FVM) [28] and finite element method (FEM) [29].In Section2, there is a detailed discussion of the model formulation and choice of accounting for the expressions.The rest of the paper is arranged as follows. Section2 explains the mathematical model. Then, for further investigation, we nondimensionalize it. Subsequently, a detailed discussion of a numerical scheme for the resulting equation is discussed in Section 3. Furthermore, a contribution is made to the Section 4, where the simulation results are presented and the effect on tumour cell population density is discussed during an interaction with CD8+ cytotoxic T-cell in its Subsection 4.6. In Section 5, we conclude our key findings and give some open problems.


## 2. Mathematical Formulation

Initially, Kuznetsov et al. made a seminal contribution to modelling tumour-immune interactions using a prey-predator system [14]. Their model studied the interactions between the tumour cell and one type of immune cell. Furthermore, El-Gohary advanced the model by studying the interaction between tumour cells and two types of immune cells and assuming that tumour cell mortality is directly proportional to tumour cell density [15]. More recently, Agrawal et al. [17] modified the El-Gohary [15] modelling approach and incorporated Holling’s type-II functional response to explain the interaction between tumour and immune cells. According to Holling’s type-II response, the mortality rate drops as tumour cell density grows, and the maximal number of tumour cells dies at low density. The model is governed by a system of ordinary differential equations with an independent time variable t. The dependent variables are the population density of the tumour cell, hunting cell, and resting cell, denoted by T, H, and R, respectively. The hunting and resting cells are part of the immune cell population termed effector cells. The system of equations is expressed as follows:(3)dTdt=q+r1T1−Tk1−α1THK1+T,(4)dHdt=α2HRK2+R−d1H,(5)dRdt=r2R1−Rk2−α2HRK2+R−d2R,where q denotes the rate at which normal cells become tumour cells. The term r1T1−T/k1 represents tumour cell proliferation, with r1 denoting its growth rate, and k1 implying its maximum carrying capacity. The term α1TH/K1+T describes the interaction of tumour cells with hunting cells, with α1 referring to the rate at which tumour cells are attacked by hunting cells and K1 referring to the half-saturation population density of tumour cells. The expression α2HR/K2+R describes the conversion of resting cells into hunting cells, with α2 implying the rate of conversion and the constant K2 implying a half-saturation population density of resting cells. The programmed cell death rate of hunting cells is denoted by d1. The term r2R1−R/k2 expresses the proliferation of resting cells, with r2 denoting the resting cell growth rate and k2 denoting the resting cell maximum carrying capacity. The programmed cell death of resting cells is denoted by d2.We extend the work of El-Gohary [15] and Agrawal et al. [17] by a few modifications discussed in the following lines. Firstly, we add diffusion to each cell’s population density; hence, the system of ordinary differential equations becomes the system of partial differential equations. Secondly, we use the BD functional response for tumour-CD8+ cytotoxic T cell interaction instead of Holling’s type-II functional response. The BD functional response makes the attacking rate of CD8+ cytotoxic T cells independent of tumour cell density, even if tumour cell density is high and describes tumour cell behaviour as dependent on an abundance of immune cells. Simultaneously, CD4+ helping T cells are supposed to promote the proliferation of CD8+ cytotoxic T cells that is expressed by Holling’s type-III functional response. Thirdly, nutrient-driven growth and death of tumour cells and the chemotaxis (directed migration towards the gradient of nutrient concentration) are incorporated in tumour cells. In addition, the activation and exhaustion of CD8+ cytotoxic T cells are expressed using Holling’s type-III functional response. Moreover, only CD8+ cytotoxic T cells are involved in killing tumour cells, whereas CD4+ helping T cells promote CD8+ cytotoxic T cell proliferation.The dependent variables for which we need to find the solution are population density of tumour cellsT, CD8+ cytotoxic T cells H, and CD4+ helping T cells R. In addition, a nutrient concentration denoted by cc>0 diffuses from the tumour’s periphery. The variable c affects the proliferation and migratory behaviour of tumour cells. The independent variables are space x and time t. Hence, we obtain the following model equation, which is a nonlinear system of the partial differential equations are given by(6)∂T∂t=DT∂2T∂x2−δ∂∂xT∂logc∂x−α1Hλ+T+γHT+qT,(7)∂H∂t=DH∂2H∂x2+α2R2h12+R2H+α3T2h22+T2H−α4T2h32+T2H−μHH,(8)∂R∂t=DR∂2R∂x2+rRR1−Rk2−μRR,with initial condition at t=0, T>0, H>0, R>0 and subject to Neumann’s boundary conditions,(9)∂T∂x=∂H∂x=∂R∂x=0atx=0andt>0,∂T∂x=∂H∂x=∂R∂x=0atx=Landt>0.In equation (6) on the right side, the first term represents random the motion of tumour cells (diffusion) with DT as the diffusion coefficient. The second term represents chemotaxis with δ as the chemotaxis coefficient. This term is chosen analogous to bacterial chemotaxis [30–33] and obeys Webber–Fechner law. The third term represents BD functional response described in equation (2), representing the interaction between tumour cells and CD8+ cytotoxic T cells with α1 denoting the maximum attacking rate of CD8+ cytotoxic T cells at tumour cells, λ as saturation constant, and γ denoting the interference parameter of CD8+ cytotoxic T cells. The term qT is accounted from [34] which is expressed as(10)qT=T1−T/k1S0c1+S1c−S2+S3c1+S4cT.The first term on the right side of equation (10) represents the tumour cell birth rate, and the second term describes the rate of tumour cell death, where S0, S1, S2, S3, and S4 are the constants that regulate the rate at which cells reproduce and die. Tumour cell birth and death are nutrient-dependent; as c approaches zero, the birth rate becomes zero. In addition, as c rises, so does the cell birth rate. When c approaches infinity in the limiting case of the first term, the birth rate becomes T1−T/k1S0c/1+S1c. Furthermore, the condition S3/S4≤S2 is imposed, ensuring that the death rate remains a decreasing function of nutrient concentration. If c=0, then the death rate term is TS2, and as c approaches infinity, the death rate term becomes TS3/S4.In equation (7), the first term on the right side is the random motion of CD8+ cytotoxic T cells with DH as the diffusion coefficient. The second term denotes the proliferation of CD8+ cytotoxic T cells due to CD4+ helper T cells with α2 as the maximum proliferation rate and h1 as the saturation constant. The third term denotes the recruitment of CD8+ cytotoxic T cells due to the presence of tumour cells with α3 as the recruitment rate and h2 as the saturation constant. The fourth term denotes the death of CD8+ cytotoxic T cells exhaustion while killing tumour cells with α4 as exhaustion rate and h3 as saturation constant. The fifth term denotes the natural death of CD8+ cytotoxic T cells with μH as the death coefficient.In equation (8), the first term is the random motion of CD4+ helper T cells with DR as diffusion coefficient. The second term represents the logistic growth of CD4+ helper T cells with rR as the growth coefficient and k2 as carrying capacity. The third term denotes the death of CD4+ helper T cells naturally with μR as the death coefficient.Nutrients required for cell growth disseminate through an external blood vessel in a tumour environment. The diffusion equation for nutrient concentration is represented as(11)∂2c∂x2=Qc,with cx,0>0 and subject to boundary conditions given by(12)∂c0,t∂x=0,cL,t=cout.The nutrient consumption rate accounted from [34] is governed by the equation given as(13)Qc=Q0Tc1+Q1c,where Q0 and Q1 are non-negative constants.


### 2.1. Nondimensionalization

On scaling, the independent and dependent variables to get corresponding dimensionless quantities denoted with a “∗” symbol are represented in the following manner:(14)t=t∗1+S1coutS0cout,x=Lx∗,T=k1T∗,H=k1H∗,R=k2R∗,c=coutc∗.Furthermore, we substitute these scaled variables in equations (6)–(8), (11), and remove the “∗” symbol to avoid notational inconvenience. The transformed equations are written as follows:(15)∂T∂t=∂2T∂x2−δ1∂∂xT∂logc∂x+T1−T1+s1c1+s1c−s2+s3c1+s4cT−a1THλ1+T+γH,(16)∂H∂t=D1∂2H∂x2+a2R2λ22+R2H+a3T2λ32+T2H−a4T2λ42+T2H−μ1H,(17)∂R∂t=D2∂2R∂x2+rR1−R−μ2R,(18)∂2c∂x2=Q0^Tc1+Q1^c.The scaled parameters derived from the nondimensionalization process are given as(19)L=DT1+S1coutS0cout,δ1=δDT,D1=DHDT,D2=DRDT,Q0^=k1L2,s1=S1cout,s2=S21+S1coutS0cout,s3=S21+S1coutS0,s4=S4cout,Q1^=Q1cout,a1=α1S0cout1+S1cout,a2=α2S0cout1+S1cout,a3=α3S0cout1+S1cout,a4=α4S0cout1+S1cout,μ1=μHS0cout1+S1cout,μ2=μRS0cout1+S1cout,r=rRS0cout1+S1cout,λ1=λk1,λ2=h12k22,λ3=h22k12,λ4=h32k12,with initial condition given as T>0, H>0, R>0, and c>0, and the two set of boundary conditions, firstly Neumann’s boundary condition,(20)∂T∂x=∂H∂x=∂R∂x=0atx=0andt>0,∂T∂x=∂H∂x=∂R∂x=0atx=Landt>0.Also, the boundary conditions for nutrient concentration are given as(21)∂c0,t∂x=0,cL,t=cout.


## 2.1. Nondimensionalization

On scaling, the independent and dependent variables to get corresponding dimensionless quantities denoted with a “∗” symbol are represented in the following manner:(14)t=t∗1+S1coutS0cout,x=Lx∗,T=k1T∗,H=k1H∗,R=k2R∗,c=coutc∗.Furthermore, we substitute these scaled variables in equations (6)–(8), (11), and remove the “∗” symbol to avoid notational inconvenience. The transformed equations are written as follows:(15)∂T∂t=∂2T∂x2−δ1∂∂xT∂logc∂x+T1−T1+s1c1+s1c−s2+s3c1+s4cT−a1THλ1+T+γH,(16)∂H∂t=D1∂2H∂x2+a2R2λ22+R2H+a3T2λ32+T2H−a4T2λ42+T2H−μ1H,(17)∂R∂t=D2∂2R∂x2+rR1−R−μ2R,(18)∂2c∂x2=Q0^Tc1+Q1^c.The scaled parameters derived from the nondimensionalization process are given as(19)L=DT1+S1coutS0cout,δ1=δDT,D1=DHDT,D2=DRDT,Q0^=k1L2,s1=S1cout,s2=S21+S1coutS0cout,s3=S21+S1coutS0,s4=S4cout,Q1^=Q1cout,a1=α1S0cout1+S1cout,a2=α2S0cout1+S1cout,a3=α3S0cout1+S1cout,a4=α4S0cout1+S1cout,μ1=μHS0cout1+S1cout,μ2=μRS0cout1+S1cout,r=rRS0cout1+S1cout,λ1=λk1,λ2=h12k22,λ3=h22k12,λ4=h32k12,with initial condition given as T>0, H>0, R>0, and c>0, and the two set of boundary conditions, firstly Neumann’s boundary condition,(20)∂T∂x=∂H∂x=∂R∂x=0atx=0andt>0,∂T∂x=∂H∂x=∂R∂x=0atx=Landt>0.Also, the boundary conditions for nutrient concentration are given as(21)∂c0,t∂x=0,cL,t=cout.


## 3. Numerical Discretization

The population density of tumour cells, CD8+ cytotoxic T cells, and CD4+ helper T cells equation is approximated using FVM, and the nutrient concentration equation is approximated using FEM. In FVM, we approximate flux using an upwind scheme. In FEM, the weak formulation for nutrient concentration equations uses continuous piecewise linear polynomials.


### 3.1. Finite Volume Scheme for Cell Phase Volume Fraction

We rearrange the population density of tumour cells, CD8+ cytotoxic T cell, and CD4+ helper T cells given in equations (15)–(21). The equation (15) is rewritten as(22)∂T∂t+vT∂T∂x=∂2T∂x2+ST,where ST and vT are given as(23)ST=T1+s11−Tc1+s1c−s2+s3c1+s4c−δ11c∂2c∂x2−1c∂c∂x2−a1Hλ1+T+γ1H,vT=δ11c∂c∂x.Also, (16) and (17) are given as(24)∂H∂t=D1∂2H∂x2+SH,(25)∂R∂t=D2∂2R∂x2+SR,where(26)SH=a2R2λ22+R2H+a3T2λ32+T2H−a4T2λ42+T2H−μ1H,SR=rR1−R−μ2R.Assuming the following initial and boundary conditions:(27)T=T0,H=H0,R=R0,att=0,∂T∂x=0,∂H∂x=0,∂R∂x=0,atx=0,∂T∂x=0,∂H∂x=0,∂R∂x=0,atx=L.Let0=x0<⋯<xM=1 and 0=t0<⋯<TN=T be the uniform discretization of domain I=0,1 and time t=0,T, respectively. We set step size ∆x=xi+1−xi and time step ∆t=tj+1−tj. We represent the ith cell of a control volume by Ii=xi−1/2,xi+1/2, where xi+1/2=xi+xi+1/2. Equations (22)–(24) is integrated in ith cell over [xi−1/2, xi+1/2 ] ×tj,tj+1, and then the backward scheme is used for time derivative. The resulting integral form yields(28)∫IiTx,tj−Tx,tj−1dx+vTi∫tj−1tjTxi+1/2−Txi−1/2dt=∫tj−1tj∂T∂xxi+1/2−∂T∂xxi−1/2dt+∫tj−1tjS¯T∆xdt,∫IiHx,tj−Hx,tj−1dx=∫tj−1tjD1∂H∂xxi+1/2−∂H∂xxi−1/2dt+∫tj−1tjS¯H∆xdt,∫IiRx,tj−Hx,tj−1dx=∫tj−1tjD2∂R∂xxi+1/2−∂R∂xxi−1/2dt+∫tj−1tjS¯H∆xdt.The right-side derivative term is approximated by(29)∂A∂xxi+1/2=Axi+1,t−Axi,t∆x,∂A∂xxi−1/2=Axi,t−Axi−1,t∆x,where A∈T,H,R. Furthermore, approximating the flux using an upwind scheme, the given discrete form of an equation is(30)Tij−1=Tij+∆t∆xvTiTij−Ti−1j−∆t∆x2Ti+1j−2Tij+Ti−1j+∆tSTi,ifvTi>0,Tij+∆t∆xvTiTi+1j−Tij−∆t∆x2Ti+1j−2Tij+Ti−1j+∆tSTi,ifvTi<0,Hij−1=Hij−D1∆t∆x2Hi+1j−2Hij+Hi−1j+∆tSHi,Rij−1=Rij−D2∆t∆x2Ri+1j−2Rij+Ri−1j+∆tSRi,where(31)Aij∼1∆x∫IiTx,tjdx.The values ofSTi, SHi, SRi and vTi are provided as follows:(32)STi=Tij−1+s11−Tijcij1+s1cij+s2+s3cij1+s4cij+δ1cijci+1j−2cij+ci−1j∆x2+1cijci+1j−ci−1j2∆x2+a1Hijλ1+Tij+γ1Hij,SHi=Hij−a2Rij2λ22+Rij2−a3Tij2λ32+Tij2+a4Tij2λ42+Tij2+μ1,SRi=Rij−r1−Rij+μ2,vTi=δ11cijci+1j−ci−1j2∆x.The above discrete system can calculate the unknown variablesT,H, and R.


### 3.2. Finite Element Method for Nutrient Concentration

Here, we discuss the variational formulations and finite element discretization of the nutrient concentration in equation (19). The computational domain is defined as I=0,1 and in given space W=v∈H1I:v1=0, where H1I is separable Hilbert space. We multiply equation (19) by a test function w∈W and apply integration by parts and given boundary conditions to obtain(33)∫I∂c∂x∂w∂xdx=−∫IwQ0^Tc1+Q1^cdx+∫∂Iw∂c∂x.Equation (19) has a weak solution is defined by c∈H10,1 such that c1=1 and equation (35) is satisfied for every w∈W. We define each element by Ii=xi−1,xi, 1≤i≤M, and P1Ii as a set of all continuous linear polynomials on Ii. The finite-dimensional subspace is given by(34)Wh≔vh:vh∈CI,vhIi∈P1Ii,1≤i≤M,vh1=0,where CI represents the space of continuous real-valued functions defined on I. The discrete problem seeks a solution of ch∈Wh with ch1=1 such that for every wh∈Wh, the unknown function c is approximated to produce the discrete form represented as(35)c=∑h=1Mchwh.For each weight functionwh, the FEM determines the residual’s weak form. The weighted residual for the given problem is(36)∫I∑h=1M∂ch∂x∂wh∂xdx+∫Iwh∑h=1MchQ0^T1+Q1^chdx−∫∂Iwh∂ch∂x=0.Although the above integral is global,wh is nonzero over the two elements that involve one node. Therefore, the local integral is expressed as(37)∑h=1M∫I∂ch∂x∂wh∂ξ+whchQ0^T1+Q1^chdx=∫∂Iwh∂ch∂x.Therefore, the matrix representation of equation (39) is(38)K+Q0^Ti1+Q1^chMch=F.The algebraic system in equation (38) can be solved to calculate ch, whereas the expressions are given below to evaluate K, M, and F for each element.(39)K=∫xh−1xh+1∂vh∂x∂wh∂xdx=1Δx1−1−11,M=∫xh−1xh+1vhwhdx=Δx62112,F=∫∂Iwh∂ch∂x=00,sinceatboundary∂ch∂x∂I=0.Here,F is a column vector.


## 3.1. Finite Volume Scheme for Cell Phase Volume Fraction

We rearrange the population density of tumour cells, CD8+ cytotoxic T cell, and CD4+ helper T cells given in equations (15)–(21). The equation (15) is rewritten as(22)∂T∂t+vT∂T∂x=∂2T∂x2+ST,where ST and vT are given as(23)ST=T1+s11−Tc1+s1c−s2+s3c1+s4c−δ11c∂2c∂x2−1c∂c∂x2−a1Hλ1+T+γ1H,vT=δ11c∂c∂x.Also, (16) and (17) are given as(24)∂H∂t=D1∂2H∂x2+SH,(25)∂R∂t=D2∂2R∂x2+SR,where(26)SH=a2R2λ22+R2H+a3T2λ32+T2H−a4T2λ42+T2H−μ1H,SR=rR1−R−μ2R.Assuming the following initial and boundary conditions:(27)T=T0,H=H0,R=R0,att=0,∂T∂x=0,∂H∂x=0,∂R∂x=0,atx=0,∂T∂x=0,∂H∂x=0,∂R∂x=0,atx=L.Let0=x0<⋯<xM=1 and 0=t0<⋯<TN=T be the uniform discretization of domain I=0,1 and time t=0,T, respectively. We set step size ∆x=xi+1−xi and time step ∆t=tj+1−tj. We represent the ith cell of a control volume by Ii=xi−1/2,xi+1/2, where xi+1/2=xi+xi+1/2. Equations (22)–(24) is integrated in ith cell over [xi−1/2, xi+1/2 ] ×tj,tj+1, and then the backward scheme is used for time derivative. The resulting integral form yields(28)∫IiTx,tj−Tx,tj−1dx+vTi∫tj−1tjTxi+1/2−Txi−1/2dt=∫tj−1tj∂T∂xxi+1/2−∂T∂xxi−1/2dt+∫tj−1tjS¯T∆xdt,∫IiHx,tj−Hx,tj−1dx=∫tj−1tjD1∂H∂xxi+1/2−∂H∂xxi−1/2dt+∫tj−1tjS¯H∆xdt,∫IiRx,tj−Hx,tj−1dx=∫tj−1tjD2∂R∂xxi+1/2−∂R∂xxi−1/2dt+∫tj−1tjS¯H∆xdt.The right-side derivative term is approximated by(29)∂A∂xxi+1/2=Axi+1,t−Axi,t∆x,∂A∂xxi−1/2=Axi,t−Axi−1,t∆x,where A∈T,H,R. Furthermore, approximating the flux using an upwind scheme, the given discrete form of an equation is(30)Tij−1=Tij+∆t∆xvTiTij−Ti−1j−∆t∆x2Ti+1j−2Tij+Ti−1j+∆tSTi,ifvTi>0,Tij+∆t∆xvTiTi+1j−Tij−∆t∆x2Ti+1j−2Tij+Ti−1j+∆tSTi,ifvTi<0,Hij−1=Hij−D1∆t∆x2Hi+1j−2Hij+Hi−1j+∆tSHi,Rij−1=Rij−D2∆t∆x2Ri+1j−2Rij+Ri−1j+∆tSRi,where(31)Aij∼1∆x∫IiTx,tjdx.The values ofSTi, SHi, SRi and vTi are provided as follows:(32)STi=Tij−1+s11−Tijcij1+s1cij+s2+s3cij1+s4cij+δ1cijci+1j−2cij+ci−1j∆x2+1cijci+1j−ci−1j2∆x2+a1Hijλ1+Tij+γ1Hij,SHi=Hij−a2Rij2λ22+Rij2−a3Tij2λ32+Tij2+a4Tij2λ42+Tij2+μ1,SRi=Rij−r1−Rij+μ2,vTi=δ11cijci+1j−ci−1j2∆x.The above discrete system can calculate the unknown variablesT,H, and R.


## 3.2. Finite Element Method for Nutrient Concentration

Here, we discuss the variational formulations and finite element discretization of the nutrient concentration in equation (19). The computational domain is defined as I=0,1 and in given space W=v∈H1I:v1=0, where H1I is separable Hilbert space. We multiply equation (19) by a test function w∈W and apply integration by parts and given boundary conditions to obtain(33)∫I∂c∂x∂w∂xdx=−∫IwQ0^Tc1+Q1^cdx+∫∂Iw∂c∂x.Equation (19) has a weak solution is defined by c∈H10,1 such that c1=1 and equation (35) is satisfied for every w∈W. We define each element by Ii=xi−1,xi, 1≤i≤M, and P1Ii as a set of all continuous linear polynomials on Ii. The finite-dimensional subspace is given by(34)Wh≔vh:vh∈CI,vhIi∈P1Ii,1≤i≤M,vh1=0,where CI represents the space of continuous real-valued functions defined on I. The discrete problem seeks a solution of ch∈Wh with ch1=1 such that for every wh∈Wh, the unknown function c is approximated to produce the discrete form represented as(35)c=∑h=1Mchwh.For each weight functionwh, the FEM determines the residual’s weak form. The weighted residual for the given problem is(36)∫I∑h=1M∂ch∂x∂wh∂xdx+∫Iwh∑h=1MchQ0^T1+Q1^chdx−∫∂Iwh∂ch∂x=0.Although the above integral is global,wh is nonzero over the two elements that involve one node. Therefore, the local integral is expressed as(37)∑h=1M∫I∂ch∂x∂wh∂ξ+whchQ0^T1+Q1^chdx=∫∂Iwh∂ch∂x.Therefore, the matrix representation of equation (39) is(38)K+Q0^Ti1+Q1^chMch=F.The algebraic system in equation (38) can be solved to calculate ch, whereas the expressions are given below to evaluate K, M, and F for each element.(39)K=∫xh−1xh+1∂vh∂x∂wh∂xdx=1Δx1−1−11,M=∫xh−1xh+1vhwhdx=Δx62112,F=∫∂Iwh∂ch∂x=00,sinceatboundary∂ch∂x∂I=0.Here,F is a column vector.


## 4. Numerical Results

This section aims to present the simulation results of a nonlinear system of partial differential equations that depicts immune cell behaviour in the tumour microenvironment in the absence of treatment. Table1 lists the parameter values used in the simulation. Some parameters marked as “—” are selected through performing numerical simulation several times. The initial population density of tumour cells, CD8+ cytotoxic T cells, and CD4+ helper T cells in the dimensionless unit taken in a column vector form is given as T0,x=0.8,…,0.8′, H0,x=0.4,…,0.4′ and R0,x=0.5,…,0.5′. Final time is T=200 and ∆t=10−3. Spatial domain x = [0, 1] is divided into 100 equal parts with step size h=10−2. The outcomes of numerical simulation are presented in furthersections.Table 1
The parameter values in dimensionless form.
ParameterValuesReferencea12.25—a30.1—D11.5—λ10.2—λ30.2—Q0^0.5[34]s1, s410[34]γ0.5—μ10.05[17]a20.05—a40.6—D21.10—λ20.1—λ40.5—Q1^0[34]S2,S30.5[34]r0.5[17]μ20.1—


### 4.1. Diffusion Effect on the Density of Tumour Cells and CD8+ Cytotoxic T Cells

The diffusion effect illustrates the nonhomogeneous spread of the population density of tumour cells and CD8+ cytotoxic T cells in a given domain [0, 1]. In the case of tumour cells, it is caused when they are overcrowded and begin to repel each other, which helps in predicting the spread of tumour cells. Consequently, it is crucial to examine how they are spread, which will aid in the treatment process. According to Figure2(a), it is observed that as DT increases, the distribution of the density of tumour cells in the region decreases. Each curve has an initial climb and a subsequent drop. Whereas, in Figure 2(b), as D1 increases, the distribution of density of CD8+ cytotoxic T cells in the region increases. Each curve starts out decreasing and then starts to increase. This indicates as tumour cells diffuse in a region, the CD8+ cytotoxic T cells target the population of tumour cells that are more densely populated and ignore the population of tumour cells with a lower density, causing the density of tumour cells to rise initially and then fall over time. While attacking the tumour cell population, these CD8+ cytotoxic T cells may increase their population density and diffuse in a region because they are activated by dying tumour cells.Figure 2
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of diffusion at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying diffusion coefficient DT in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying diffusion coefficient D1 in a given space x. The rest of the parameters are taken from Table 1.
(a)(b)Now, suppose that the density and spatial distribution influence the diffusion of tumour cells. Here, we describe the diffusion coefficient as a nonlinear space-dependent and density-dependent function:(40)DT=r1TTm+1x2,where r1 and m are the positive constants. To ensure that the diffusion is slow when the cell density is either too low or too high, we choose m>1.Figure3, graph shows that the tumour cell density diffusion decreases over time. This is because, according to our assumption, diffusion depends on density and spatial distribution. As time increases, density of tumour cell decreases due to other governing factors like limited nutrients or an attack by CD8+ cytotoxic T cells. As a result, the diffusion of tumour cells also reduces in a region.Figure 3
An example of tumour cell density due to density-dependent and space-dependent diffusion coefficientDTT,x at various time intervals, i.e., t=0.1, t=0.5 and t=1. Parameter values are r1 = 1.8 and m = 4.


### 4.2. Effect of Maximum Proliferation Parametera2 of CD8+ Cytotoxic T Cells due to CD4+ Helper T Cells

Here, the effect of the maximum proliferation parametera2 can be observed in Figure 4. It indirectly impacts the population density of tumour cells. As CD4+ helper T cells also promote the growth of CD8+ cytotoxic T cells, the parameter a2 is directly propositional to the population density of CD8+ cytotoxic T cells. Thus, the population density of CD8+ cytotoxic T cells increases with an increase in a2 as shown in Figure 4(b). Whereas an increase in population density of CD8+ cytotoxic T cells impacts the density of tumour cells inversely. As a result, the population density of tumour cells decreases with an increase in a2 as shown in Figure 4(a).Figure 4
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of maximum proliferation parametera2 at time t=1. Each curve in (a) shows an evolution of the density of tumour cells on varying maximum proliferation parameter a2 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying maximum proliferation parameter a2 in a given space x.
(a)(b)


### 4.3. Effect of Maximum Recruitment Parametera3 of CD8+ Cytotoxic T Cells due to Dying Tumour Cells

Further, the effect of the maximum recruitment parametera3 can be observed in Figure 5. In Figure 5(a), as a3 increases, the distribution of tumour cell density decreases in the region. Whereas, in Figure 5(b), as a3 increases, the distribution of CD8+ cytotoxic T cells density increases in the region. By comparing 5(a) and 5(b), it is observed that when the density of tumour cells is low, the density of CD8+ cytotoxic T cells is high, and vice versa, suggesting that the death of tumour cells enhances the recruitment of CD8+ cytotoxic T cells and that recruitment is low at high tumour cell density.Figure 5
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of maximum recruitment parameter at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying maximum recruitment parameter a3 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying maximum recruitment parameter a3 in a given space x.
(a)(b)


### 4.4. Effect of Exhaustion Parametera4 of CD8+ Cytotoxic T Cells due to Attacking on Tumour Cells

In addition, the effect of the exhaustion parametera4 can be observed in Figure 6. In Figure 6(a), as a4 increases, the distribution of tumour cell density increases in the region. Whereas, in Figure 6(b), as a4 increases, the distribution of CD8+ cytotoxic T cells density decreases in the region. This indicates that when tumour cell density is low, CD8+ cytotoxic T cells are less exhausted and retain a high density; however, as tumour cell density grows, CD8+ cytotoxic T cell exhaustion increases, resulting in a drop in density as observed in Figures 6(a) and 6(b).Figure 6
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of maximum exhaustion parameter at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying maximum exhaustion parameter a4 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying maximum exhaustion parameter a4 in a given space x.
(a)(b)


### 4.5. Effect of Saturation Constantsλ3 and λ4

Furthermore, the effect of the saturation constantλ3 can be observed in Figure 7. In Figure 7(a), as λ3 increases, the distribution of tumour cell density increases in the region which indicates that there are not many tumour cells dying. As a result, CD8+ cytotoxic T cells do not activate in response to dying tumour cells and their density declines as λ3 increases, as shown in 7(b).Figure 7
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of saturation constant at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying saturation constant λ3 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying saturation constant λ3 in a given space x.
(a)(b)Similarly, the effect of the saturation constantλ4 can be observed in Figure 8. In Figure 8(a), as λ4 increases, the distribution of tumour cell density increases in the region, signifying that fewer cancer cells are dying. As a result, CD8+ cytotoxic T cells do not get exhausted, thus preserving their own population density. Therefore, as seen in Figure 8(b), as λ4 increases, the distribution of CD8+ cytotoxic T cells density increases in the region. Here, the population density of both cells increases, but the density of tumour cells increases more rapidly with an increase in time.Figure 8
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of saturation constant at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying saturation constant λ4 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying saturation constant λ4 in a given space x.
(a)(b)Additionally, more simulations are performed to compare the population density of each cell.Finding the exact solution of the model equation is an intricate task as the equations are coupled and highly nonlinear. To demonstrate the effectiveness of the applied numerical scheme, we will compute the error norms∙1, ∙2, and ∙∞. The error norms are defined as follows:(41)∙1=∑j=1NThj−Th∗j,∙2=h∑j=1NThj−Th∗j2,∙∞=maxjThj−Th∗j,where Th is the approximated solution with step size h and Th∗ is the numerical solution with step size h∗=2h. We calculate the convergence ratio, which is defined as the ratio of successive errors of the given norms for the proposed scheme. If the ratio approaches the numbers 2, 4, 6, and so on, it indicates that the order of convergence is 1, 2, and 3. The details about this approach can be found in [35].Figure 9
The graph shows the population density of tumour cells, CD8+ cytotoxic T cells and CD4+ helper T cells with respect to spacex and for the parameters given in Table 1 at time t=1.From Table2, as the ratio is approaching 2, the computed solution has first-order accuracy.Table 2
Error norms for the solution of T given in Figure9 at t=1.
MethodN.1r.2r.∞rFVM-FEM50.0016—0.0020—0.0034—108.8500e − 041.800.00111.810.00191.7204.5900e − 041.925.5652e − 041.970.00101.9402.3300e − 041.962.8296e − 041.965.2000e − 041.92801.1662e − 041.991.4172e − 041.992.6000e − 0421605.8813e − 051.987.1567e − 051.981.4000e − 041.85Figure 10
For the parameters given in Table1 at time t=200, the graph shows (a) space vs. population density and (b) time vs. population density plot.
(a)(b)From Table3, as the ratio is approaching 2, the computed solution has first-order accuracy.Table 3
Error norms for the solution of T given in Figure10 at t = 200.
MethodN.1r.2r.∞rFVM-FEM51.3178e − 38—1.5994e − 38—2.5980e − 38—106.2003e − 392.127.3331e − 392.181.2320e − 382.10203.0182e − 392.053.5238e − 392.086.0000e − 392.05401.4899e − 392.021.7284e − 392.032.9600e − 392.02807.4020e − 402.018.5603e − 402.011.4700e − 392.011603.6896e − 402.004.2610e − 402.007.3000e − 402.01All the graphical results are computed in the first order of approximation in time and second order approximation in space.


### 4.6. Discussion

The outcomes of numerical simulation are discussed in this section. The parameter values are chosen so that the population density of tumour cells can be minimized or eliminated. Besides, in equations (15)–(17), the rates of activation a3 and saturation constants λ3, as well as the rates of exhaustion of a4 and saturation constants λ4, are chosen in such a way that activation occurs at lower antigen levels than exhaustion, and the latter dominates at high antigen levels. As a result, a3<a4 and λ3<λ4. While performing the simulation, we observed that the parameter values a2 and a3 should be minimum because increasing these values for a long-time period causes overactivation of CD8+ cytotoxic T cells, leading to autoimmune diseases such as arthritis and diabetes. This can be observed in Figures 5 and 6, where at time t = 1, the population density of CD8+ cytotoxic T cells increases rapidly with an increase in a2 and a3.Moreover, from Figure2, we can say that the diffusion of tumour cells can be controlled if we increase the diffusion of CD8+ cytotoxic T cells in a limited manner. As shown in Figure 9, the population density of each cell at time t=1 coexists in space. The population density of tumour cells is lower than its initial population density at time t=1. This is because CD8+ cytotoxic T cells attack tumour cells. When CD8+ cytotoxic T cells attack tumour cells, they get exhausted, resulting in a population drop from its original density. Overall, CD8+ cytotoxic T cell population density is lower than tumour cell population density. This is because the tumour cells are acquiring enough nutrients to survive at t=1; therefore, their birth rate is higher than mortality due to CD8+ cytotoxic T cell attack or nutrient deficiency. However, as time increases, the population density of tumour cells decreases to zero in a given space, as seen in Figure 11(a). Whereas after a certain time t=15, the population density of CD8+ cytotoxic T cells and CD4+ helper T cells reaches a saturation threshold as observed in Figures 11(b) and 11(c). The long-time behaviour of each cell population density at time t=200 is given in Figure 10, indicating that using the parameter values in Table 1, the tumour cell density can be eliminated.Figure 11
The simulation result for the parameters given in Table1 till the time t=20. Graph (a) represents the population density of tumour cells, (b) represents the population density of CD8+ cytotoxic T cells and (c) represents the population density of CD4+ helper T cells, with respect to space x in varying time t.
(a)(b)(c)In addition, as shown in Figures12(a) and 12(b), an increase in the exhaustion rate of CD8+ cytotoxic T cells causes the elimination of CD8+ cytotoxic T cells in a given space and time, as well as an increase in the population density of cancer cells. The population density of tumours decreases at first, then increases to a certain point before saturating due to limited nutrients, as illustrated in the time vs. population density graph Figure 12(b). Biologically, as antigen levels rise, CD8+ cytotoxic T cells are activated first, followed by CD8+ cytotoxic T cell exhaustion. In the next section, we conclude the results and discussion.Figure 12
The graph shows (a) space vs. population density and (b) time vs. population density plot when the value ofa4=0.8 at time t=200.
(a)(b)


## 4.1. Diffusion Effect on the Density of Tumour Cells and CD8+ Cytotoxic T Cells

The diffusion effect illustrates the nonhomogeneous spread of the population density of tumour cells and CD8+ cytotoxic T cells in a given domain [0, 1]. In the case of tumour cells, it is caused when they are overcrowded and begin to repel each other, which helps in predicting the spread of tumour cells. Consequently, it is crucial to examine how they are spread, which will aid in the treatment process. According to Figure2(a), it is observed that as DT increases, the distribution of the density of tumour cells in the region decreases. Each curve has an initial climb and a subsequent drop. Whereas, in Figure 2(b), as D1 increases, the distribution of density of CD8+ cytotoxic T cells in the region increases. Each curve starts out decreasing and then starts to increase. This indicates as tumour cells diffuse in a region, the CD8+ cytotoxic T cells target the population of tumour cells that are more densely populated and ignore the population of tumour cells with a lower density, causing the density of tumour cells to rise initially and then fall over time. While attacking the tumour cell population, these CD8+ cytotoxic T cells may increase their population density and diffuse in a region because they are activated by dying tumour cells.Figure 2
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of diffusion at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying diffusion coefficient DT in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying diffusion coefficient D1 in a given space x. The rest of the parameters are taken from Table 1.
(a)(b)Now, suppose that the density and spatial distribution influence the diffusion of tumour cells. Here, we describe the diffusion coefficient as a nonlinear space-dependent and density-dependent function:(40)DT=r1TTm+1x2,where r1 and m are the positive constants. To ensure that the diffusion is slow when the cell density is either too low or too high, we choose m>1.Figure3, graph shows that the tumour cell density diffusion decreases over time. This is because, according to our assumption, diffusion depends on density and spatial distribution. As time increases, density of tumour cell decreases due to other governing factors like limited nutrients or an attack by CD8+ cytotoxic T cells. As a result, the diffusion of tumour cells also reduces in a region.Figure 3
An example of tumour cell density due to density-dependent and space-dependent diffusion coefficientDTT,x at various time intervals, i.e., t=0.1, t=0.5 and t=1. Parameter values are r1 = 1.8 and m = 4.


## 4.2. Effect of Maximum Proliferation Parametera2 of CD8+ Cytotoxic T Cells due to CD4+ Helper T Cells

Here, the effect of the maximum proliferation parametera2 can be observed in Figure 4. It indirectly impacts the population density of tumour cells. As CD4+ helper T cells also promote the growth of CD8+ cytotoxic T cells, the parameter a2 is directly propositional to the population density of CD8+ cytotoxic T cells. Thus, the population density of CD8+ cytotoxic T cells increases with an increase in a2 as shown in Figure 4(b). Whereas an increase in population density of CD8+ cytotoxic T cells impacts the density of tumour cells inversely. As a result, the population density of tumour cells decreases with an increase in a2 as shown in Figure 4(a).Figure 4
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of maximum proliferation parametera2 at time t=1. Each curve in (a) shows an evolution of the density of tumour cells on varying maximum proliferation parameter a2 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying maximum proliferation parameter a2 in a given space x.
(a)(b)


## 4.3. Effect of Maximum Recruitment Parametera3 of CD8+ Cytotoxic T Cells due to Dying Tumour Cells

Further, the effect of the maximum recruitment parametera3 can be observed in Figure 5. In Figure 5(a), as a3 increases, the distribution of tumour cell density decreases in the region. Whereas, in Figure 5(b), as a3 increases, the distribution of CD8+ cytotoxic T cells density increases in the region. By comparing 5(a) and 5(b), it is observed that when the density of tumour cells is low, the density of CD8+ cytotoxic T cells is high, and vice versa, suggesting that the death of tumour cells enhances the recruitment of CD8+ cytotoxic T cells and that recruitment is low at high tumour cell density.Figure 5
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of maximum recruitment parameter at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying maximum recruitment parameter a3 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying maximum recruitment parameter a3 in a given space x.
(a)(b)


## 4.4. Effect of Exhaustion Parametera4 of CD8+ Cytotoxic T Cells due to Attacking on Tumour Cells

In addition, the effect of the exhaustion parametera4 can be observed in Figure 6. In Figure 6(a), as a4 increases, the distribution of tumour cell density increases in the region. Whereas, in Figure 6(b), as a4 increases, the distribution of CD8+ cytotoxic T cells density decreases in the region. This indicates that when tumour cell density is low, CD8+ cytotoxic T cells are less exhausted and retain a high density; however, as tumour cell density grows, CD8+ cytotoxic T cell exhaustion increases, resulting in a drop in density as observed in Figures 6(a) and 6(b).Figure 6
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of maximum exhaustion parameter at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying maximum exhaustion parameter a4 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying maximum exhaustion parameter a4 in a given space x.
(a)(b)


## 4.5. Effect of Saturation Constantsλ3 and λ4

Furthermore, the effect of the saturation constantλ3 can be observed in Figure 7. In Figure 7(a), as λ3 increases, the distribution of tumour cell density increases in the region which indicates that there are not many tumour cells dying. As a result, CD8+ cytotoxic T cells do not activate in response to dying tumour cells and their density declines as λ3 increases, as shown in 7(b).Figure 7
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of saturation constant at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying saturation constant λ3 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying saturation constant λ3 in a given space x.
(a)(b)Similarly, the effect of the saturation constantλ4 can be observed in Figure 8. In Figure 8(a), as λ4 increases, the distribution of tumour cell density increases in the region, signifying that fewer cancer cells are dying. As a result, CD8+ cytotoxic T cells do not get exhausted, thus preserving their own population density. Therefore, as seen in Figure 8(b), as λ4 increases, the distribution of CD8+ cytotoxic T cells density increases in the region. Here, the population density of both cells increases, but the density of tumour cells increases more rapidly with an increase in time.Figure 8
Plots show the density of tumour cells and CD8+ cytotoxic T cells spread due to the effect of saturation constant at timet=1. Each curve in (a) shows an evolution of the density of tumour cells on varying saturation constant λ4 in a given space x, and each curve in (b) exhibits an evolution of the density of CD8+ cytotoxic T cells on varying saturation constant λ4 in a given space x.
(a)(b)Additionally, more simulations are performed to compare the population density of each cell.Finding the exact solution of the model equation is an intricate task as the equations are coupled and highly nonlinear. To demonstrate the effectiveness of the applied numerical scheme, we will compute the error norms∙1, ∙2, and ∙∞. The error norms are defined as follows:(41)∙1=∑j=1NThj−Th∗j,∙2=h∑j=1NThj−Th∗j2,∙∞=maxjThj−Th∗j,where Th is the approximated solution with step size h and Th∗ is the numerical solution with step size h∗=2h. We calculate the convergence ratio, which is defined as the ratio of successive errors of the given norms for the proposed scheme. If the ratio approaches the numbers 2, 4, 6, and so on, it indicates that the order of convergence is 1, 2, and 3. The details about this approach can be found in [35].Figure 9
The graph shows the population density of tumour cells, CD8+ cytotoxic T cells and CD4+ helper T cells with respect to spacex and for the parameters given in Table 1 at time t=1.From Table2, as the ratio is approaching 2, the computed solution has first-order accuracy.Table 2
Error norms for the solution of T given in Figure9 at t=1.
MethodN.1r.2r.∞rFVM-FEM50.0016—0.0020—0.0034—108.8500e − 041.800.00111.810.00191.7204.5900e − 041.925.5652e − 041.970.00101.9402.3300e − 041.962.8296e − 041.965.2000e − 041.92801.1662e − 041.991.4172e − 041.992.6000e − 0421605.8813e − 051.987.1567e − 051.981.4000e − 041.85Figure 10
For the parameters given in Table1 at time t=200, the graph shows (a) space vs. population density and (b) time vs. population density plot.
(a)(b)From Table3, as the ratio is approaching 2, the computed solution has first-order accuracy.Table 3
Error norms for the solution of T given in Figure10 at t = 200.
MethodN.1r.2r.∞rFVM-FEM51.3178e − 38—1.5994e − 38—2.5980e − 38—106.2003e − 392.127.3331e − 392.181.2320e − 382.10203.0182e − 392.053.5238e − 392.086.0000e − 392.05401.4899e − 392.021.7284e − 392.032.9600e − 392.02807.4020e − 402.018.5603e − 402.011.4700e − 392.011603.6896e − 402.004.2610e − 402.007.3000e − 402.01All the graphical results are computed in the first order of approximation in time and second order approximation in space.


## 4.6. Discussion

The outcomes of numerical simulation are discussed in this section. The parameter values are chosen so that the population density of tumour cells can be minimized or eliminated. Besides, in equations (15)–(17), the rates of activation a3 and saturation constants λ3, as well as the rates of exhaustion of a4 and saturation constants λ4, are chosen in such a way that activation occurs at lower antigen levels than exhaustion, and the latter dominates at high antigen levels. As a result, a3<a4 and λ3<λ4. While performing the simulation, we observed that the parameter values a2 and a3 should be minimum because increasing these values for a long-time period causes overactivation of CD8+ cytotoxic T cells, leading to autoimmune diseases such as arthritis and diabetes. This can be observed in Figures 5 and 6, where at time t = 1, the population density of CD8+ cytotoxic T cells increases rapidly with an increase in a2 and a3.Moreover, from Figure2, we can say that the diffusion of tumour cells can be controlled if we increase the diffusion of CD8+ cytotoxic T cells in a limited manner. As shown in Figure 9, the population density of each cell at time t=1 coexists in space. The population density of tumour cells is lower than its initial population density at time t=1. This is because CD8+ cytotoxic T cells attack tumour cells. When CD8+ cytotoxic T cells attack tumour cells, they get exhausted, resulting in a population drop from its original density. Overall, CD8+ cytotoxic T cell population density is lower than tumour cell population density. This is because the tumour cells are acquiring enough nutrients to survive at t=1; therefore, their birth rate is higher than mortality due to CD8+ cytotoxic T cell attack or nutrient deficiency. However, as time increases, the population density of tumour cells decreases to zero in a given space, as seen in Figure 11(a). Whereas after a certain time t=15, the population density of CD8+ cytotoxic T cells and CD4+ helper T cells reaches a saturation threshold as observed in Figures 11(b) and 11(c). The long-time behaviour of each cell population density at time t=200 is given in Figure 10, indicating that using the parameter values in Table 1, the tumour cell density can be eliminated.Figure 11
The simulation result for the parameters given in Table1 till the time t=20. Graph (a) represents the population density of tumour cells, (b) represents the population density of CD8+ cytotoxic T cells and (c) represents the population density of CD4+ helper T cells, with respect to space x in varying time t.
(a)(b)(c)In addition, as shown in Figures12(a) and 12(b), an increase in the exhaustion rate of CD8+ cytotoxic T cells causes the elimination of CD8+ cytotoxic T cells in a given space and time, as well as an increase in the population density of cancer cells. The population density of tumours decreases at first, then increases to a certain point before saturating due to limited nutrients, as illustrated in the time vs. population density graph Figure 12(b). Biologically, as antigen levels rise, CD8+ cytotoxic T cells are activated first, followed by CD8+ cytotoxic T cell exhaustion. In the next section, we conclude the results and discussion.Figure 12
The graph shows (a) space vs. population density and (b) time vs. population density plot when the value ofa4=0.8 at time t=200.
(a)(b)


## 5. Conclusion

In conclusion, our model effectively predicts tumour cell eradication at avascular stages of tumour cell population density with immune cell populations. The findings show that decreasing the strength of CD8+ cytotoxic T cells (i.e. increasing the exhaustion rate of CD8+ cytotoxic T cells) causes tumour cells to escape from immune cells. Furthermore, enhancing the interaction strength of T cells (i.e. the activation rate of CD8+ cytotoxic T cells) to attack tumour cells could significantly reduce tumour population density and even lead to the elimination of tumour cells. Moreover, varying diffusion parameters can slow down a tumour’s growth. Thus, the activation rate, exhaustion rate of CD8+ cytotoxic T cells, and diffusion coefficients (DT, D1) are significant parameters for studying the spatial distribution of the population density of tumour cells.Moreover, understanding the spatial dynamics of immune and tumour cells is crucial for developing effective treatments. As mentioned in the introduction, some mathematical models analyze the spatial and temporal dynamics of the interaction between the immune system and the tumour by including more complex phenomena like the concentration of chemokines or adding more variables. Still, the chemotaxis of tumour cells was neglected. We incorporated this phenomenon, which helps determine the spread of tumour cells due to available nutrients. If the spread is more rapid, any treatment to treat the tumour will take longer to effect. In our model, the immune cells attacking tumour cells are expressed by BD functional response, which describes the behaviour of tumour cells is dependent on an abundance of immune cells. In contrast, earlier models used expression in which the behaviour of tumour cells is independent of immune cells when immune cells were attacking tumour cells. The present model in equations (15)–(17) is discretized using FVM and equation (19) for nutrient concentration is discretized using FEM. The obtained discrete equation is solved using MATLAB version R2020a. We have also calculated error norms to get the order of accuracy. Because we used a first-order numerical scheme, the findings demonstrated that the current technique provides first-order precision for modelling the tumour-immune interaction.The tumour-immune system interaction cannot be predicted further when the immune cell population density is exhausted. This exhaustion occurs due to an attack on tumour cells, and this can limit further predictions. In order to avoid such circumstances, an optimized dosage of immunotherapy or a combination of therapies at regular intervals is needed that will maintain the population density of CD8+ cytotoxic T cells and CD4+ helper T cells. The optimized dosage is vital because the overactivation of CD8+ cytotoxic T cells can cause autoimmune diseases such as arthritis and diabetes. In future work, we can extend the model by incorporating the optimized drug dosage. We can also consider the role of other types of immune systems, like the innate immune system, which will provide more insight into the mechanism involved in tumour-immune interaction. Furthermore, because the immune system plays an essential role in removing viruses and other pathogens, the present model helps understand other biological phenomena, such as wound healing problems and investigating the current coronavirus outbreak 2019 (COVID-19).



---

*Source: 9006678-2023-01-03.xml*