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.
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.
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 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 aswhere denotes the half-saturation constant, is the immune cell’s maximal attacking rate, and 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 aswhere 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 Section 2, 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. Section 2 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 . The dependent variables are the population density of the tumour cell, hunting cell, and resting cell, denoted by , , and , respectively. The hunting and resting cells are part of the immune cell population termed effector cells. The system of equations is expressed as follows:where denotes the rate at which normal cells become tumour cells. The term represents tumour cell proliferation, with denoting its growth rate, and implying its maximum carrying capacity. The term describes the interaction of tumour cells with hunting cells, with referring to the rate at which tumour cells are attacked by hunting cells and referring to the half-saturation population density of tumour cells. The expression describes the conversion of resting cells into hunting cells, with implying the rate of conversion and the constant implying a half-saturation population density of resting cells. The programmed cell death rate of hunting cells is denoted by . The term expresses the proliferation of resting cells, with denoting the resting cell growth rate and denoting the resting cell maximum carrying capacity. The programmed cell death of resting cells is denoted by .
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 cells , CD8+ cytotoxic T cells , and CD4+ helping T cells . In addition, a nutrient concentration denoted by c diffuses from the tumour’s periphery. The variable c affects the proliferation and migratory behaviour of tumour cells. The independent variables are space and time . Hence, we obtain the following model equation, which is a nonlinear system of the partial differential equations are given bywith initial condition at , , , and subject to Neumann’s boundary conditions,
In equation (6) on the right side, the first term represents random the motion of tumour cells (diffusion) with 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 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 is accounted from [34] which is expressed as
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 , , , , and are the constants that regulate the rate at which cells reproduce and die. Tumour cell birth and death are nutrient-dependent; as approaches zero, the birth rate becomes zero. In addition, as rises, so does the cell birth rate. When approaches infinity in the limiting case of the first term, the birth rate becomes . Furthermore, the condition is imposed, ensuring that the death rate remains a decreasing function of nutrient concentration. If , then the death rate term is , and as approaches infinity, the death rate term becomes .
In equation (7), the first term on the right side is the random motion of CD8+ cytotoxic T cells with as the diffusion coefficient. The second term denotes the proliferation of CD8+ cytotoxic T cells due to CD4+ helper T cells with as the maximum proliferation rate and as the saturation constant. The third term denotes the recruitment of CD8+ cytotoxic T cells due to the presence of tumour cells with as the recruitment rate and as the saturation constant. The fourth term denotes the death of CD8+ cytotoxic T cells exhaustion while killing tumour cells with as exhaustion rate and as saturation constant. The fifth term denotes the natural death of CD8+ cytotoxic T cells with as the death coefficient.
In equation (8), the first term is the random motion of CD4+ helper T cells with as diffusion coefficient. The second term represents the logistic growth of CD4+ helper T cells with as the growth coefficient and as carrying capacity. The third term denotes the death of CD4+ helper T cells naturally with 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 aswith and subject to boundary conditions given by
The nutrient consumption rate accounted from [34] is governed by the equation given aswhere and 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:
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:
The scaled parameters derived from the nondimensionalization process are given aswith initial condition given as , , , and , and the two set of boundary conditions, firstly Neumann’s boundary condition,
Also, the boundary conditions for nutrient concentration are given as
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 aswhere and are given as
Also, (16) and (17) are given aswhere
Assuming the following initial and boundary conditions:
Let and be the uniform discretization of domain and time , respectively. We set step size and time step . We represent the cell of a control volume by , where . Equations (22)–(24) is integrated in cell over [, ] , and then the backward scheme is used for time derivative. The resulting integral form yields
The right-side derivative term is approximated bywhere . Furthermore, approximating the flux using an upwind scheme, the given discrete form of an equation iswhere
The values of , , and are provided as follows:
The above discrete system can calculate the unknown variables , and .
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 and in given space , where is separable Hilbert space. We multiply equation (19) by a test function and apply integration by parts and given boundary conditions to obtain
Equation (19) has a weak solution is defined by such that and equation (35) is satisfied for every . We define each element by , , and as a set of all continuous linear polynomials on . The finite-dimensional subspace is given bywhere represents the space of continuous real-valued functions defined on . The discrete problem seeks a solution of with such that for every , the unknown function is approximated to produce the discrete form represented as
For each weight function , the FEM determines the residual’s weak form. The weighted residual for the given problem is
Although the above integral is global, is nonzero over the two elements that involve one node. Therefore, the local integral is expressed as
Therefore, the matrix representation of equation (39) is
The algebraic system in equation (38) can be solved to calculate , whereas the expressions are given below to evaluate , , and for each element.
Here, 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. Table 1 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 , and . Final time is and . Spatial domain = [0, 1] is divided into 100 equal parts with step size . The outcomes of numerical simulation are presented in furthersections.
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 Figure 2(a), it is observed that as 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 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.
(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:where and are the positive constants. To ensure that the diffusion is slow when the cell density is either too low or too high, we choose .
Figure 3, 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.
4.2. Effect of Maximum Proliferation Parameter of CD8+ Cytotoxic T Cells due to CD4+ Helper T Cells
Here, the effect of the maximum proliferation parameter 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 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 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 as shown in Figure 4(a).
(a)
(b)
4.3. Effect of Maximum Recruitment Parameter of CD8+ Cytotoxic T Cells due to Dying Tumour Cells
Further, the effect of the maximum recruitment parameter can be observed in Figure 5. In Figure 5(a), as increases, the distribution of tumour cell density decreases in the region. Whereas, in Figure 5(b), as 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.
(a)
(b)
4.4. Effect of Exhaustion Parameter of CD8+ Cytotoxic T Cells due to Attacking on Tumour Cells
In addition, the effect of the exhaustion parameter can be observed in Figure 6. In Figure 6(a), as increases, the distribution of tumour cell density increases in the region. Whereas, in Figure 6(b), as 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).
(a)
(b)
4.5. Effect of Saturation Constants and
Furthermore, the effect of the saturation constant can be observed in Figure 7. In Figure 7(a), as 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 increases, as shown in 7(b).
(a)
(b)
Similarly, the effect of the saturation constant can be observed in Figure 8. In Figure 8(a), as 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 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.
(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 , , and . The error norms are defined as follows:where is the approximated solution with step size and is the numerical solution with step size . 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].
From Table 2, as the ratio is approaching 2, the computed solution has first-order accuracy.
(a)
(b)
From Table 3, as the ratio is approaching 2, the computed solution has first-order accuracy.
All 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 and saturation constants , as well as the rates of exhaustion of and saturation constants , 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, and . While performing the simulation, we observed that the parameter values and 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 and .
Moreover, from Figure 2, 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 coexists in space. The population density of tumour cells is lower than its initial population density at time . 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 ; 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 , 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 is given in Figure 10, indicating that using the parameter values in Table 1, the tumour cell density can be eliminated.
(a)
(b)
(c)
In addition, as shown in Figures 12(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.
(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 (, ) 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).
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.