Reliable A Posteriori Error Estimator for a Multi-scale Cancer Invasion Model

A cell is known to be the fundamental unit of life in all living organisms. For the organisms to operate accurately, these cells must develop and split in a controlled manner according to a specific set of rules. Generally, cells are immobile and can develop, multiply, and kill themselves in a self-regulated mode. Occasionally, cells may expand freely without considering the standard ratio of growth and death. Cells damage host tissue if they develop and reproduce by forgetting the body’s needs and constraints. In a broad sense, cancer is a condition in which normal cells begin to replicate uncontrollably. It is a category of disease that has numerous origins and evolves over time and space. Therefore, cancer is a complex phenomenon to foresee. Comprehending the mechanism of cancer progression is critical to detecting and treating this disease. Understanding the process of cancer evolution relies on numerical simulations to a great extent, which entrusts us to see the growth and spreading of the tumor. Cancer models are treated as reaction-diffusion equations in numerous works of literature [6].

There is ample literature for a posterior analysis of finite element methods for elliptic and parabolic equations, to note a few [22, 3, 12]. The robust, reliable, and efficient residual-based error estimate for the singularly perturbed reaction-diffusion problem was derived in [21, 4, 2, 20]. In [22], reliable and efficient residual-based error estimators are discussed for linear and nonlinear elliptic equations and linear and quasilinear parabolic equations. In [12], the duality approach is used for a posteriori error estimation for a system of reaction-diffusion equations. To the author’s knowledge, theoretical a-posterior error estimates for the cancer invasion model were not done in the past. Also, previous works focused less on a coupled system of nonlinear parabolic equations to derive the theoretical estimates, where one unknown depends on the gradient of the other. The basic idea behind our work can be traced back to the method of residual-based error estimation in [22].

Numerical techniques to solve the tumor invasion model and associated continuum mathematical models are discussed in the following. Solving a mathematical model of angiogenesis by using the finite difference scheme was presented in [5, 10]. In [15], the research delved into the examination of the four-species tumor growth model through the application of a two-dimensional mixed Finite Element Method (FEM). Meanwhile, [23] presented numerical findings for the continuum tumor invasion model, employing the FEM while considering the expansion of capillaries within a two-dimensional spatial domain. To enhance computational efficiency in practical real-world scenarios, Adaptive Mesh Refinement (AMR) techniques, whether in the temporal or spatial domain, have found favorable utilization, as discussed in references [11, 19, 18, 9, 8]. In the study documented in [16], various time discretization methods and the adaptive finite volume approach were employed for conducting two-dimensional simulations. Furthermore, [24] harnessed an adaptive Finite Element Method (FEM) to address equations related to tumor angiogenesis in a two-dimensional context. Applied adaptive FEM to investigate a simplified three-equation reaction model that computes aspects of tumor-induced angiogenesis deterministically in [18]. The study presented numerical results for geometric models across dimensions, specifically considering $d-$dimensional geometry, where $d=\{1,2,3\}$ and utilized explicit time-stepping schemes. Recently, an administration technique for the bookkeeping of AMR on (hyper-)rectangular meshes is presented in [17] and a computational framework of space-time adaptivity for this mode in our previous work [7]. In [7], a computational framework based on the parallel space-adaptive techniques was presented based on the gradient type error estimator for spatial discretization. Still, the rigorous theoretical estimates of the error estimator were not presented. The previous studies investigated cancer or related models, and none of the works showed the theoretical study of a-posteriori error estimates and their numerical converge studies of cancer invasion models. This is the main focus of the current paper.

Let $u(x,t)$, $v(x,t)$, and $w(x,t)$ be the three unknowns describing the cancer cell’s density, extracellular matrix (ECM) density, and concentration of matrix-degrading enzymes (MDE), respectively. The governing equations for the dimensionless form of the cancer invasion model are as follows.

where $d_{1}(u,v,w)$ is the density-dependent nonlinear diffusion coefficient, ${\chi_{u}(v)}$ is the haptotactic sensitivity function, $d_{2}$ is diffusion constant, and $\lambda,\rho,\eta,\alpha,\beta$ are positive constants.
Assuming the interactions of the cancer cells with ECM and the degradation of ECM by MDE take place in an isolated system, zero-flux type boundary conditions are imposed on the model.

where $\xi$ is the outward normal vector on $\partial\Omega$. The initial conditions for the given system of equations (1) are

with $\epsilon=0.01$ and $r^{2}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}$.

In this current work, we focus on the theoretical framework of a posteriori error estimates using the residual-based error estimator, and their numerical realization is investigated. This is one of the main goals of this paper. In section 2, the weak formulation of the model problem (1) and its space and temporal discretization are discussed. In section 3, first, we derive a relationship between the error and residuals for the given problem. And then, a residual-based error estimate is proposed by finding an upper bound for the residuals. The numerical results are discussed in section 4.

Let $\Omega\subset\mathbb{R}^{3}$ be the Lipschitz continuous bounded domain. As usual, we denote $H^{1}(\Omega)$ as the Sobolev space of the functions in $L^{2}(\Omega)$ whose weak derivatives also belong to $L^{2}(\Omega)$. To define spaces involving time, let $X$ be a real Banach space equipped with norm $\|\cdot\|$ and $T(>0)\in\mathbb{R}$. The space $L^{p}(0,T;X)$ consists of all measurable functions $u~{}:~{}[0,T]\rightarrow X$ with

The space $H^{1}(0,T;X)$ consists of all functions $u\in L^{2}(0,T;X)$ such that $u^{\prime}$ exists in the weak sense and belongs to $L^{2}(0,T;X)$.

Let $X(\Omega)$ be a Banach space and $X^{*}$ be its dual. Let $\|\cdot\|_{*}$ denote the dual norm of elements in $X^{*}$. So, from the definition of the dual norm, we have

where $<\cdot,\cdot>$ denotes the duality pairing between $X$ and $X^{*}$. We refer [13, 1] for more details on Sobolev spaces and Bochner spaces. To derive the weak formulation of model problem (1), we multiply the equations (1) by test functions and integrate over $\Omega$. After applying integration by parts, we get the weak formulation of the model problem:

For given $u_{0},v_{0},w_{0}\in L^{2}(\Omega)$ find $u,v,w\in L^{2}(0,T;H^{1}(\Omega))$ with $\displaystyle\frac{\partial u}{\partial t},\frac{\partial v}{\partial t},\frac{\partial w}{\partial t}\in L^{2}(0,T;(H^{1}(\Omega))^{*})$ such that

for all $\psi_{1},\psi_{2},\psi_{3}\in H^{1}(\Omega)$. Further, $(H^{1}(\Omega))^{*}$ is the dual space of $H^{1}(\Omega))$.

To define the space discretization of weak formulation (2), we adopt the conforming Galerkin method. Let $\Omega_{h}$ be the shape regular triangulation of $\Omega$ and $V_{h}$ be finite-dimensional subspace of $H^{1}(\Omega)$, where $h$ is the discretization parameter. The discrete problem reads as : find $u_{h},v_{h},w_{h}\in V_{h}$ such that for a.e. $t\in(0,T)$ and for all $\psi_{1,h},\psi_{2,h},\psi_{3,h}\in V_{h}$

For each $(t_{i},t_{j})\subset(0,T)$, consider the function space

We apply the backward Euler method for time discretization. Let $0=t^{0}<t^{1}<t^{2}<\cdots<t^{J}=T$ be a partition of the considered time interval $[0,T]$ with step size $\tau^{n}=t^{n}-t^{n-1},~{}n=1,2,\ldots,J$. Thus it follows from (2)

where $u^{n+1}_{h}(t),v^{n+1}_{h}(t),w^{n+1}_{h}(t)\in V_{h}$ for $t\in[n\tau^{n},(n+1)\tau^{n}),1\leq n\leq J-1$. To derive the relationship between the error and residuals, we have the following assumptions on nonlinear diffusion and sensitivity function of $\eqref{Model:dim}$.

For all $s_{1},s_{2},s_{3},\tilde{s}_{1},\tilde{s}_{2},\tilde{s}_{3}\in\mathbb{R}$,

For the existence and uniqueness of the fully discretized system (2), we refer to [7].
Here, we provide the finite element setting to derive the error estimates. To approximate the model problem (1), we associate a shape regular triangulation $(\mathcal{T}_{h}^{n})_{h>0}$ of $\Omega$ with each time step $t_{n}$. We denote the set of all faces of the partition $(\mathcal{T}_{h}^{n})$ as $\mathcal{E}_{h}^{n}$. Furthermore, let $\tilde{\mathcal{T}}_{h}^{n}$ be the refinement of both triangulation $\mathcal{T}_{h}^{n}$ and $\mathcal{T}_{h}^{n+1}$. Also, we denote by $h_{K}$ the diameter of an element $K\in\tilde{\mathcal{T}}_{h}^{n}$ and $h_{K^{\prime}}$ denotes the diameter of an element $K^{\prime}\in\mathcal{T}_{h}^{n}$ such that $K\subseteq K^{\prime}$.

This section is devoted to the study of residual operators and error estimation for the discretized problem (2). In subsection 3.1, we define the residual operators and an upper bound for the error is derived in terms of the given initial data and the residuals. In subsection 3.2, we derive the upper bound for the residuals by splitting them into spatial and temporal residuals.