Stability analysis and simulations of tumor growth model based on system of reaction-diffusion equation in two-dimensions

1.   Introduction

Partial differential equations (PDEs) play a crucial role in mathematical modeling across various scientific and engineering disciplines. For example, Maxwell's equations are used to govern how electric and magnetic fields evolve. The Navier-Stokes equations describe the motion of fluid substances and are essential in aerodynamics, weather modeling, and oceanography. The heat equation models the distribution and evolution of temperature in a given region over time. PDEs help in understanding stress and strain distribution within materials, crucial for building bridges, buildings, and various structures. Reaction-diffusion equations model the spread of a species in an environment and interactions between different species. Further, these PDEs are highly effective in simulating the dynamics of tumor progression. A prevalent method involves utilizing a reaction-diffusion equation, which delineates the alteration in the concentration of certain substances due to local reactions (like chemical reactions, growth, and decay) coupled with diffusion (the dispersion stemming from random movement). Within modeling tumor growth, the focus is often on the cell density of tumor cells as the substance of interest. Here, the reaction component accounts for the proliferation rate of tumor cells (such as through mitosis), while the diffusion component reflects the movement and invasion of tumor cells into surrounding tissues ^[1,2].

The mechanism behind carcinogenic mutations remains intricate and to date is not yet fully understood, despite the extensive focus of many researchers on specific signaling pathways associated with it. For instance, there is a wealth of studies examining the role of genetic mutations. Furthermore, compared to other processes related to tumor growth, the mathematical modeling of this particular process is not as developed. Recent times have witnessed significant research interest in the mathematical modeling of brain tumor growth, especially concerning reaction-diffusion models. The utilization of these models has proven crucial in gaining a more profound understanding of brain tumor growth patterns, which has furthered the advancement of tailored treatment approaches ^{[3,4,5,6,7,8,9,10]}.

The ongoing quest to find an effective and lasting treatment for tumors remains a significant challenge for scientists. However, considerable advancements have been made in developing new methodologies that have proven successful in reducing or even eliminating tumors. Among these, mathematical modeling stands out as a crucial tool for enhancing cancer therapy. Carcinogenesis is a highly complex process that requires extensive study for a comprehensive understanding. Tumors originate from one or more normal cells that have experienced malignant transformation. The immune system's response to tumors varies based on the tumor's antigenicity. A cell with substantial mutations forms a tumor more easily recognized as foreign (or more antigenic) compared to one that only slightly deviates from a healthy cell ^[11]. For various cancer types, the carcinogenic process can be divided into different stages, typically ranging from $4$ to $7$, depending on the tumor type^[12]. In ^[13], the author examines a system of Lotka-Volterra type delay differential equations (DDEs) with diffusion that models the mutation of cells from normal to malignant. This model is primarily based on mutations in three different environmental conditions: favorable, competitive, and unfavorable. They concentrated on identifying running wave solutions. Subsequently, a broader analysis was conducted on models with two and an arbitrary number of mutation steps in ^[14]. The case of unfavorable conditions often seen as a consequence of treatments like chemotherapy. Instead of detailing each mutation stage, the analysis simplifies the process by incorporating a suitable time delay, effectively transforming the system of $n + 1$ ordinary differential equations with diffusion into a pair of equations with time delay. A comparative study is undertaken between the dynamics of the ODE system, the system with delay, and the system with both delay and diffusion. Analytically, it is demonstrated that the stability of the positive steady state (when present) is highly dependent on the delay's magnitude. Furthermore, researchers have explored the potential for stability switches in steady states and have demonstrated that a change in the stability of the positive steady state is improbable in the absence of delay in the model.

In ^[15], the authors introduced spatial heterogeneity into their diffusion model by assigning different diffusion rates, utilizing a reaction-diffusion equation to depict the tumor cell density. While the exact ratio of diffusion coefficients between these brain regions can vary, they specifically postulated a five fold difference in diffusion rates. A detailed representation of white and gray matter within a synthetic brain for MR simulations was constructed, allowing the observation and analysis of unique growth patterns arising from the varied diffusion rates in the brain.

Subsequent progress was reported in ^[16], where the model was enhanced to include the anisotropic expansion of gliomas, inspired by the tendency of glial cells to migrate along fiber tracts. Comparing their simulation results with two clinical cases, it was demonstrated that the anisotropic diffusion model more accurately represented the shape and growth kinetics of low-grade gliomas near the insula compared to isotropic diffusion models. This suggests the importance of considering cell diffusion tensor anisotropy in simulating glioma growth and spread. Additionally, the study in ^[17] developed a mathematical model based on diffusion tensor imaging to identify major white matter tracts in the brain. The anisotropic model presented slight superiority over other reaction-diffusion models, particularly for tumors exhibiting strong anisotropy. A new approach for estimating parameters in reaction-diffusion tumor growth models using medical images taken at different times was introduced in ^[18]. This technique computes patient-specific parameters of the model, enabling the adjustment of the general model to fit individual patient data. The research indicated that, when the tumor cell proliferation rate was held constant, several other parameters could be uniquely determined. This methodological advancement offers a more personalized and accurate representation of tumor growth, enhancing the potential for tailored patient treatment strategies.

Numerical solutions of the reaction-diffusion equation model for tumor growth are essential for developing novel cancer therapies, strengthening patient-specific treatment plans, and expanding our knowledge of tumor dynamics. This strategy offers a potent tool for addressing the difficulties associated with cancer research and treatment by fusing theoretical understanding with real-world applications. The interrelationship of multiple biological processes, such as cell migration, proliferation, and interaction with their surroundings, characterizes tumors as extremely complex systems. These dynamics are modeled by the reaction-diffusion equation, which simulates the development of cancer cells (reaction) and their movement through the tissue (diffusion) along with other chemicals and nutrients. This makes it possible to analyze tumor growth patterns and how the environment affects the tumor's development in great detail. The reaction-diffusion model's numerical solutions allow for the prediction of a tumor's growth and possible spread of the tumor to other sections of the body. Recognizing these trends can aid in forecasting the course of the illness and in creating more potent treatment plans to stop the spread of cancer cells. Treatment scenarios, including as chemotherapy, radiation therapy, and innovative therapeutic techniques, can be simulated using numerical models. Before implementing these treatments in clinical settings, researchers can find possibly viable medications and optimize treatment regimens by evaluating how these treatments affect tumor growth inside the model.

The primary objective of this work is to investigate theoretically as well as numerically a two-stage tumour growth model in a two-dimensional based on system of reaction-diffusion equations by employing the Euler method in time and the Fourier spectral method in space with appropriate initial and boundary conditions. Detailed stability analysis conditions are derived using the fixed points of the system. The choice of using the Fourier spectral method is driven by their high accuracy and minimal phase error in solving the specified problem. A significant advantage of spectral methods is their exponentially decaying error, implying an infinite convergence rate in space even with a minimal number of collocation points using the fast Fourier transform. These methods offer flexibility in choosing basis functions, which can be selected based on the problem at hand. However, their implementation can be challenging, mainly because they use global functions as basis functions. This characteristic makes them less suitable for handling local features and sharp gradients, as seen in phenomena, like Gibbs phenomenon ^{[19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38]}. Due to their ability to describe the intricate dynamics of tumor growth and dissemination, differential equations and their fractional form are essential tools in the understanding and treatment of fatal diseases such as cancer. They support the modeling of interactions between cancer cells and their surroundings, the forecasting of disease development, and the assessment of the efficacy of different approaches to treatment. To improve patient outcomes, targeted medicines and personalized medicine approaches are designed with the support of this mathematical framework, which offers crucial insights into the mechanisms underlying the development of cancer^{[39,40,41,42,43,44,45,46,47,48,49]}.

The structure of this paper is as follows: Section 2 introduces the governing equations of the reaction-diffusion system with detailed implementation of the Fourier method in space and the Euler method in time. Section 3 evaluates the stability of the model equation using the equilibrium points, followed by numerical discussions in Section 4. Finally, Section 5 provides a concise conclusion of the findings.

2.   Formulation of the mathematical model and description of the numerical method

Mathematical models, especially those utilizing PDEs, are highly beneficial for simulating tumor progression. A prevalent technique involves reaction-diffusion equations that depict the variation in the concentration of one or several substances due to local reactions (such as chemical reactions, growth, and decay) coupled with diffusion (the distribution caused by random motion). Specifically, in tumor growth modeling, the focus is often on the cell density of tumor cells as the substance of interest. The reaction component in this context indicates the proliferation rate of tumor cells (like through cell division), while the diffusion component illustrates the dispersion of tumor cells within the tissue.

Let $Y_j$ denotes the density of cells which is mutant at $j$-th stage at some position, where $j = 0, 1, 2, ..., n$. To develop a multi-stage model having mutant cell densities $Y_j$ at stage $j = 0, 1, 1, ..., n$, with $j = 0$ represent the first stage. The intermediate stage is occurring for $j = 1, ..., n-1$ and $j = n$ gives us the final stage. The density function can be described by the system of equations ^[50]:

where $a_j$ rate of logistic growth, $K_j$ corresponds to carrying capacity, $\Delta$ is the Laplacian operator, $\mu_j$ and $\eta_j$ represent the interaction coefficients for the corresponding stages and $D_j$ is the co-efficient of diffusion. To simplify the model equation (2.1), we can rescale the model such that $K_j \leftrightarrow 1$ for $j = 1, \ldots, n-1$, while for the last mutation stage we assume $K_n \to \infty$. After rescaling, we get:

In model equation (2.1), we denote the initial (normal) stage $j = 0$ by $u_0 = Y_0$ and the final stage $Y_n = u_n$ and assuming that the $j$-th stage is not the initial (normal) or final stage and call it the intermediate or benign stage. Then, the tumor growth model which describes a process of carcinogenesis mutations involving $n$ different stages of mutations (transitioning from normal to malignant cells) is given by:

$F$ is a particular function which is contingent on the environmental conditions surrounding tumor growth and development and whose value depends on three different environments corresponding to the the final stage $j = n$:

● In favourable conditions:

● In a competitive environment:

● In unfavourable conditions:

In model equation (2.3), the natural choice for the contingent function $F$ is the unfavourable conditions.

Consider the two-stage model in unfavourable conditions with included delay and diffusion corresponds to initial, intermediate and malignant stages. The primary contribution we make is expanding the model from ^[51,52] to two-dimensional space. After rescaling, we get:

where $u_0$ represents the density of normal cells at time $t$ and position $(x, y)$, $u_1$ and $u_2$ stand for benign and malignant cells, respectively, all constants are positive, delays $\tau_i$, $i = 1$, $2$ are non-negative, and $(x, y)\in[0, \pi]$. We consider the zero flux boundary conditions:

and initial condition:

For time discretization of model equation (2.7), we use the Euler method, while for the spacial discretization, we use the Fourier spectral method.

The Fourier spectral method is a technique for solving PDE, that is particularly useful when the problem is defined on a periodic domain. It is a type of spectral method, where solutions are approximated using a sum of functions from a certain function space. In the case of the Fourier spectral method, this space is composed of trigonometric functions (sines and cosines), which are the eigenfunctions of linear operators with periodic boundary conditions. It uses sine and cosine functions as the basis for the spectral expansion. The solution to the PDE is represented as a sum of these basis functions, each multiplied by a coefficient. Unlike finite difference or finite element methods, which are local methods, the Fourier spectral method considers the global information of the problem, offering very accurate solutions with fewer basis functions for smooth problems. For smooth problems, the method has spectral accuracy, meaning the error decreases exponentially with the number of basis functions used. It is particularly powerful for problems with periodic boundary conditions and in situations where high accuracy is needed over the entire domain. Its efficiency diminishes for problems with non-periodic boundary conditions or with sharp gradients or discontinuities.

Suppose $u(x)$ is a function defined on a periodic domain $x \in [0, L]$. $u(x)$ can be represented as a Fourier series in the form:

Here, $\hat{u}_n$ are the Fourier coefficients, and $e^{i2\pi nx/L}$ are the complex exponential that form the basis of the function space. For real functions, one uses sines and cosines, but the complex exponential form is more convenient for mathematical manipulations. The Fourier coefficients $\hat{u}_n$ capture the amplitude of the function at different frequencies and are given by:

The resultant PDE is then transformed into a set of algebraic equations for the coefficients $\hat{u}_n$. Solving these equations gives the Fourier coefficients, which you can then use to reconstruct the solution. A Discrete Fourier Transform (DFT) is used for this purpose. For a set of discrete points $u_j = u(x_j)$ where $j = 0, 1, ..., N-1$ and $x_j = jL/N$, the DFT is defined as:

where, $\hat{U}_k$ are the discrete Fourier coefficients, which approximate $\hat{u}_n$. For reconstructing the original function from its Fourier coefficients, one can use the inverse DFT:

The Fast Fourier Transform (FFT) is an efficient way to compute the DFT and its inverse. The most common FFT algorithm is the Cooley-Tukey algorithm, which recursively breaks down the DFT of any composite size $N = N_1N_2$ into many smaller DFTs, exploiting symmetries and reducing the computational cost from $O(N^2)$ to $O(N\log N)$. Then the spacial domain $[0, \pi] \times [0, \pi]$ is discretized into a grid of $N \times N$ points, where $N$ is a power of 2 for efficient FFT computation. For each function $u_i(t, x, y)$, we compute the DFT to obtain the Fourier coefficients $\hat{u}_i(t, k_x, k_y)$. This is done using the FFT. The derivative of a function in Fourier space is obtained by multiplying its Fourier coefficients by $ik$, where $k$ is the wave number. For the two-dimensional case:

The wave numbers $k_x$ and $k_y$ are related to the grid points and the size of the spatial domain. For time discretization, let $t_n = n\Delta t$ where $\Delta t$ is the time step size. The Euler method updates the solution from $t_n$ to $t_{n+1}$ as:

where $F(\hat{u}_{i, n, m}^n, t_n)$ represents the right-hand side of the discretized PDEs.

The final discretization form of model equation (2.7) using the Euler scheme in time and the Fourier spectral method in space is given by:

The non-linear terms in the model equation (2.7) are handled in real space before being transformed back into Fourier space. The system includes terms that depend on the solution at earlier times ($t-\tau_1$ and $t-\tau_2$), a buffer initialization for each $u_i$ to store the necessary past states.

This discretization turns the continuous PDEs into a system of algebraic equations, which is then solved using traditional techniques for solving linear systems of equations.

3.   Stability analysis

In this section, we will outline the stability analysis for the proposed model, Eq (2.7), employing certain numerically defined parameter values that hold biological relevance, in conjunction with the determined equilibrium states. These critical points play an essential role in assessing the population dynamics of normal, benign, and tumor cells. The stability of these equilibrium points is ascertained by numerically calculating the eigenvalues of the Jacobian matrix associated with the model. An equilibrium point is deemed stable if all the eigenvalues are negative, indicating a scenario where normal, benign, and tumor cells maintain a balanced coexistence without significant growth or decline. Conversely, an equilibrium point is considered unstable if any eigenvalue is positive, suggesting a potential for tumor cells to rapidly overtake normal cells following minor fluctuations in cell densities, potentially leading to malignancy. Our numerical simulations incorporate diffusion factors, and we demonstrate that these do not affect the stability of the cell densities. A stable equilibrium ensures that any minor deviations in cell densities due to delays eventually returns to their initial steady states. In contrast, an unstable equilibrium may result in significant deviations from steady states upon slight perturbations, raising the possibility of malignancy development.

This stability analysis is confined to cases without diffusion and delay in the model, as detailed in Eq (2.7). However, the subsequent section will extend the analysis to include these factors in our numerical simulations. For a comprehensive understanding of the model's behavior with both delay and diffusion in a one-dimensional setting, we direct the reader to ^[51,52]. These findings are also applicable to our two-dimensional case.

In the absence of both delay and diffusion, Eq (2.7) gives a system of ordinary differential equations of the form:

The Jacobian matrix for system (2.2) is given by:

The characteristic equation $p(\lambda)$ of the Jacobian matrix is:

The system of model equation (2.7) has the following six equilibrium points:

$(i)\; A_1(0, 0, 0)$ - Trivial equilibrium where all populations are extinct.

$(ii)\; A_2(0, 1, 0)$ - Only the benign cells ($u_1$) are present at their carrying capacity.

$(iii)\; A_3\left(0, \frac{1}{\phi_2}, \frac{\alpha_1\phi_2 - \alpha_1}{\phi_2\psi_2}\right)$ - Only the benign and malignant cells are present, and their populations depend on the interaction parameters $\phi_2$, $\alpha_1$, and $\psi_2$.

$(iv)\; A_4(1, 0, 0)$ - Only the normal cells ($u_0$) are present at their carrying capacity.

$(v)\; A_5\left(\frac{\alpha_0\phi_2 - \psi_1}{\alpha_0\phi_2}, \frac{1}{\phi_2}, \frac{\alpha_0\alpha_1\phi_2 - \alpha_0\alpha_1 + \alpha_0\phi_1\phi_2 - \phi_1\psi_1}{\alpha_0\phi_2\psi_2}\right)$ - A mixed equilibrium where all three populations are present. The populations depend on multiple parameters of the system.

$(vi)\; A_6\left(\frac{\alpha_1(\alpha_0 - \psi_1)}{\alpha_0\alpha_1 + \phi_1\psi_1}, \frac{\alpha_0\alpha_1 + \alpha_0\phi_1}{\alpha_0\alpha_1 + \phi_1\psi_1}, 0\right)$- An equilibrium where normal and benign cells are present, and their populations depend on the interaction parameters.

These equilibrium points represent different possible steady states of the system of model equation (2.7), where the rate of change of each cell type is zero. The actual stability and feasibility (whether the populations are non-negative and make sense biologically) of these points would depend on the specific values of the model parameters $\alpha_0, \alpha_1, \psi_1, \psi_2, \phi_1$, and $\phi_2$, which we discuss in the subsequent section. After finding the Jaocobain using these equilibrium points and using some parameter values involved in the model equations, the following observations are made:

$\diamond\; A_1(0, 0, 0)$-the trivial steady states exist and do not depends on the parameters involved in the model equation.

$\diamond\; A_2(0, 1, 0)$-it reflect that benign cells are present and a half-trivial steady exists and again is independent of the parameters.

$\diamond\; A_3\left(0, \frac{1}{\phi_2}, \frac{\alpha_1\phi_2 - \alpha_1}{\phi_2\psi_2}\right)$- a positive steady state exists for $\alpha_1 = 0.5, \psi_2 = 5, \phi_2 = 3.5$, but there is one stability switch for $\tau_2 = 1.253670632, \tau_1 = 0$.

$\diamond\; A_4(1, 0, 0)$ - a half-trivial steady state exists in the presence of normal cells ($u_0$).

$\diamond\; A_5\left(\frac{\alpha_0\phi_2 - \psi_1}{\alpha_0\phi_2}, \frac{1}{\phi_2}, \frac{\alpha_0\alpha_1\phi_2 - \alpha_0\alpha_1 + \alpha_0\phi_1\phi_2 - \phi_1\psi_1}{\alpha_0\phi_2\psi_2}\right)$ - a positive steady state exists for $\alpha_0 = 0.5; \alpha_1 = 0.3; \psi_1 = 2; \psi_2 = 5; \phi_1 = 3; \phi_2 = 5;$ and one stability switch for $\tau_2 = 1.253670632; \tau_1 = 0.$

$\diamond\; A_6\left(\frac{\alpha_1(\alpha_0 - \psi_1)}{\alpha_0\alpha_1 + \phi_1\psi_1}, \frac{\alpha_0\alpha_1 + \alpha_0\phi_1}{\alpha_0\alpha_1 + \phi_1\psi_1}, 0\right)$- a half-trivial steady state exists for $\alpha_0 = 2; \alpha_1 = 0.3; \psi_1 = 1.2; \psi_2 = 5; \phi_1 = 3; \phi_2 = 0.5;$ and stability switch with respect to $\tau_1$ for $\tau_1 = 2.726991542; \tau_2 = 0.$

4.   Numerical discussions

In this section, numerical simulations using different parameter values for the suggested model proposed in Eq (2.7) are presented ^[52]. These parameter values are chosen based on the stability analysis results. In numerical simulation, we include both delay and diffusion to see their effects. These simulations further enhance our theoretical justifications. The set of parameters without diffusion and with diffusion are given in Tables 1 and 2, respectively. We start our simulations using parameters set 1, that is, there is no diffusion and delay as shown in Figure 1, and all the cells are stable. In the absence of diffusion, it has been shown that the stability changes due to one of the delay terms using parameters sets 1, 2, and 3. This effect can be seen in Figures 2–5. The positive steady state $A_6$ exists and it does not change at all when both delays are equal to zero. It has been further observed that changing the value of a parameter may switch the steady states, especially if we change the parameter $\phi_2$ as shown in Table 1, parameter set 2, and set 3. In this case, the positive steady case disappears and is replaced by the half-steady case $A_5$. These states can be observed in Figures 3 and 4. In addition to these parameter values, we use the following initial conditions for parameter sets 1–6, in addition to the zero-flux boundary conditions given in Eq (2.8):

We then switch our simulations by including diffusion in our simulations using parameter set 6 and 8. As shown in Figures 6 and 8, diffusion does not destabilize the steady states and the cause of switch of stability is again due to one of the delay. We use the following initial conditions for simulation using set 7 for Figure 7 and set 8 for Figure 8:

Similarly, we use the below initial conditions for parameter set 9:

and

The purpose of choosing these initial conditions is to make sure that the initial cells are healthy and to see how normal cells are absorbed in the mutant and the cancer cells by a small perturbation in some parameters involved in the model equation. These effects can be seen in Figure 9. Increasing the magnitude of the delay terms may cause some small oscillation as shown in Figure 10.

5.   Conclusions

It is the aim of this work to use a robust and efficient numerical scheme for the approximate solution of tumor growth model for the justifications of theoretical claims. We studied a the model without and with delay and diffusion. It was found that the inclusion of diffusion does not effect the global steady states. Including diffusion with large coefficients may cause some oscillating behavior of the solution, but the dynamics are changed as one or both of the delay terms increases. Consequently, the intrinsic kinetic dynamics shift to an oscillatory state due to a Hopf bifurcation when the delay reaches certain values. Spatially non-homogeneous periodic solutions is obtained in the system when both delay and diffusion are present. It has been observed that adding diffusion to the model equation corresponds to global steady states, and the other parameters, particularly the delay terms are responsible for the transition in the steady states. In the future, we want to extend this methodology to the fractional version of the tumor growth model based on the reaction-diffusion system.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

This work was supported by the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia [Grant No. 5615].

Conflict of interest

The authors declare that they have no conflicts of interest.