Mathematical modeling and numerical simulation of a multiscale cancer invasion of host tissue

1.   Introduction

Mathematical modeling of biological systems especially on cancer has in recent times received much attention because, mathematical models of relevant biological processes and its associated numerical simulation, has reduced the complicated and costly experimental procedures.This is buttressed by the fact that simulation of cancer behavior across multiple biological scales in space and time is increasingly being recognized as a powerful tool that enables more accurate predictions, see ^{[1,2,3,4,6,23]} and references therein. The primary objective of multidisciplinary efforts in cancer research from the fields of mathematics, physics, and chemistry just to mention a few, is to understand the etiology of the disease and to develop suitable diagnosis and treatment.

Cancer is a generic term given to a group of diseases in which abnormal cells grow uncontrollably. Additionally, cancer cells can invade surrounding tissues and further spread to other organs in the body (http//www.who.int/health-topics). Cells become cancerous from a series of molecular (gene mutations) events that alter their fundamental properties and structural rearrangements. The resultant effects of this structural rearrangement is a characteristic abnormality in the process regulating the cells carefully regulated order of proliferation, differentiation, and death. If this abnormality continues after many years, the cells evolve through different premalignant stages to become invasive and degrade the extracellular matrix (ECM) or surrounding tissue. Such abnormal cells proliferate exponentially to form a lump called a tumor. There are two types of tumors: benign and malignant tumors, and whiles the former is not harmful the later is cancerous. Malignant tumor cells grow uncontrollably and may spread to other organs and neighboring tissues in the body through the metastasis process. Cancer development has three distinct stages: avascular, vascular, and metastatic ^[22]. A review of the literature shows that though all the three stages have been well studied and documented, work on avascular tumors dominates the number of research among the three. This does not in any way imply the importance of avascular tumors in cancer research but rather it is the simplest to model mathematically and to validate using experimental results (eg. from multicellular spheroids). Secondly, the worth of information drawn from research on avascular tumors can increase our understanding of the last two stages (vascular and metastasis). It is worth mentioning that the invasion begins as soon as the tumor cells can penetrate the surrounding tissue. Cells in the avascular tumor stage are not cancerous but form a lump due to some abnormality in the growth and agammaegation of proliferating cells.

Though the early stages of tumors are difficult to study due to their small size, education on the disease and technological advancements have contributed immensely to the study of the early stages of the tumor through in vitro experiments using multicellular spheroids, -an experimental process of growing seed cell(s) taken from a tumor cell line in a liquid medium supplemented with sufficient nutrients, see ^[13,20,21].

Tumor cells receive nutrients (eg., oxygen and glucose) supply from surrounding tissues. These nutrients reach the central core of the tumor by diffusing through layers of cells. Therefore, as the tumor increases in size, the amount of nutrients from the tumor microenvironment (TME) reaching cells at the inner rim (center) decreases from the outer layer until it becomes insufficient to sustain viable cells at the center of the tumor. Experimental studies of multicellular spheroids exhibit three distinct characteristic layers of nutrient-rich proliferating cells and an inner layer of dead (necrotic) nutrient-starved cells. These two layers are separated by a layer of nutrient deficient quiescent (hypoxic) cells. Quiescent cells are alive and can diffuse but are non-proliferating. Cells in the quiescent state may proliferate when they receive sufficient nutrients. The tumor size may be limited by the amount of nutrients that reach the inner layer of the tumor. Thus as the tumor grows in size, the volume of necrotic cells at the central core also increases. Tumor growth stops when the balance of cell proliferation and cell death is reached ^[15,22]. Therapies for cancer treatment usually target proliferating cells therefore, the population of quiescent cells may be liable to treatment failures ^[20,21,22]. It is therefore very crucial to explore and investigate the population of quiescent cells as we report in this study.

Insufficient supply of nutrients into the inner rim of the lump of cells is a contributing factor that causes the cells to initiate a process to acquire their blood supply (angiogenesis). Once the cells succeed to acquire their blood supply, the tumor cells can escape from one primary site to a secondary site to initiate secondary tumors through the circulatory system (metastasis) at other parts of the body. At the angiogenesis stage, the tumor is vascular and there is a decrease in the cell-cell and cell-matrix adhesion ^[1]. Vascular and metastatic stages are those that reduce the quality of an individual patient's life.

In furtherance to cell-cell and cell-ECM adhesion, cell surface receptors bind various ECM fibers to enable cell-matrix adhesion ^[12]. By adhering to the ECM fibres, cells can migrate from a region of lower concentration of ECM adhesive fibers to a region of higher concentration of ECM adhesive fibers in a process known as haptotaxis. Furthermore, cells migrate in response to nutrients/chemoattractant (or chemorepellent) in a process known as chemotaxis. See a detailed discussion of these processes in ^[3,6,10].

The past five years have seen tremendous in vitro and in vivo research efforts aimed at understanding the etiology of cancer cell's tissue invasion and to develop appropriate therapies for treatment. Mathematical models on cancer have either been discrete or continuum for cell dynamics ^[33]. The first mathematical model for tumors formulated in terms of densities of live and dead cells was proposed by Ward and Kind (1997, 1999), in contrast to earlier models which were formulated based on discrete cell populations. Though clinically, the latter formulation is artificial and not supported by experimental evidence. In Ward and King (1997), tumor growth was entirely due to the generation of new cell division. The analysis of interdependence between viable proliferating cell rim, growth velocities, and model parameters were also reported in their model.

In ^[11], the authors report that the relationship between solid tumors, ECM and multiple cell types in their TME enables the solid tumor to transform systematically into an invasive, metastatic and lethal disease.

Cancer cells secrete proteolytic enzymes (metalloproteinases (MMP) and or urokinase plasmogen activator (uPA)) which bind the ECM fibers to transform their conformation and configuration. Furthermore, the binding can break the natural and physical barriers that restrict the cancer cells from spreading into surrounding tissues. Mathematical models for the biologically multiscale character of cancer invasion of host tissues comes with challenges emanating from the complex interplay between the tissue-scale (macro) and the subcellular, microscale events involving the proteolytic enzymes binding to the ECM, see ^{[11,18,24,25,26,27]}. The micro and macro scales differ in space and time as we discuss in section 2. The ECM acts as a scaffold for cells as well as a signaling pathway through which the cells communicate specific intracellular and extracellular instructions for cell proliferation, migration, differentiation and most importantly, the secretion of proteolytic enzymes that destroy the ECM ^[31,32]. In their two-scale moving boundary model, the authors in ^[11] report on the effect of proteolytic enzymes (uPA) at the molecular level on the tissue-scale (macro) dynamics at the edge of the tumor invasion. Further research on cancer invasion of host tissues and investigation on effect of uPA and MMP on the evolution of cancer cells through mathematical modeling backed by experimental data include ^{[3,11,28,29,30]}. Chaplain and Lolas ^[28] used a system of reaction-diffusion-taxis (RDT) partial differential equations to explore and report on the role of uPA's in the invasion, growth and spread of cancer cells whiles Chaplain and Deakins ^[30] studied the role of MMP's in cancer growth and spread. The proof of global existence and unique classical solution for the system of RDT PDE's using a priori estimates to the model proposed by Chaplain and Lolas ^[6] is given in ^[16]. Further models for cancer invasion include ^{[1,2,3,4,5,6,7,8,9,10,17]}, and references therein.

The paper is structured as follows: In Section 2 we discuss the formulation of the mathematical model, followed in Section 3 by a nondimentionalization and scaling of the model. In Section 4, detailed discretization of the model, special treatment of the nonlinear diffusion and haptotaxis terms and numerical solutions of the model at specific time intervals is presented. In Sections 5 and 6 we discuss the findings and conclusion of the model respectively. }

2.   Formulation of the mathematical model

The study focuses on the first step of cancer invasion of ECM. We use a continuum approach to capture the entire micro and macro scale volumetric tumor growth dynamics. The model has two parts: first, a set of ordinary differential equations that models the biochemical dynamics of the tumor cells interaction with ECM at the subcellular (micro) level. In the second part, a system of partial differential equations is used to model the evolution of proliferating cells, quiescent cells, necrotic cells, and ECM densities at the macroscopic level. Unless otherwise mentioned the terms tumor cells and cancer cells are used interchangeably. The above can best be described as an advection-reaction-diffusion (ARD) system of equations. In the next subsection, we present the events at the microscopic level followed by the macroscopic level events.

2.1.   Subcellular events

Cancer cell invasion of tissues involves biochemical processes taking place at the subcellular, microscopic-scale level. These subcellular level events affect the cell motility and haptotaxis of tumor cells and intrinsically the invasion of ECM. The events at the subcellular level take place on scales in space and time that are different from the macroscale events. To observe the effects of the subcellular events at the macroscopic level, there is a time lag between the biochemical event at the subcellular level and when it is observed at the macroscopic (tissue) level. We will account for the time lag in a contractivity function used to couple the subcellular, microscopic scale events to the macroscopic scale events. It is important to mention that in between these two time-scales (micro and macro) is an intermediate time-scale known as the meso-scale (not discussed in this work), see ^[3] for more details. From the law of mass action, we can summarize the reaction between the proteolytic enzymes (l) and the ECM fibers (v) in a reversible reaction mechanism in Eq. (2.1). Equation (2.1) holds since none of the reacting proteins and products formed is volatile, precipitate or is consumed in other reactions.

where $k^+$ and $k^{-}$ are association and dissociation constants respectively. Let $l_0$ and $y$, denote the initial concentration of proteolytic enzymes and the concentration of proteolytic enzymes bound to ECM proteins ($lv$) respectively, then it holds $l_0 \in \mathcal{R}_+$. The concentration of unbound proteolytic enzymes on the cell surface is given by $l_0-y$. An ordinary differential equation (ODE) for the bound proteolytic enzymes is given in Eq. (2.2).

where $\theta$ is the microscopic time scale.

2.2.   Macroscale events

It is assumed that the onset of tumor invasion of host tissue is predominantly dominated by random motion, and haptotaxis due to spatial gradients as a result of ECM fiber degradation. At the macroscopic level, we formulate our model in terms of a continuum of densities of necrotic cells, quiescent cells, proliferating cells and ECM fibers denoted by $u(x, t)$, $q(x, t)$, $p(x, t)$ and $v(x, t)$, respectively. The model is supplemented with nutrient concentration $c(x, t)$ from the underlying tissue. The spatial coordinate in the one-dimensional domain is $x$ and $t$ is time. By this formulation, our model is geared towards an in vivo than an in vitro as this approach also allows us to account for the onset of invaded tissue in the model.

As in ^[6,34], a conservation equation for these two contributions to the tumor cell flux $J_n$ is given by $J_n = J_{random}+ J_{hapto.}$. We assume that proliferating and quiescent cells migrate whiles necrotic cells are immotile. The diffusion of cancer cells into the surrounding tissues is given by $\Phi_j(\kappa, p, q, v)\nabla j$, $j\in\{p, q\}$, where the non-linear density-dependent diffusion function $\Phi_j(\kappa, p, q, v)$ depends on the contractivity function $\kappa$, proliferating and quiescent cell densities (accounting for contact inhibition of migrating cells) and the ECM density $v(x, t)$. The cells need to adhere to adhesion sites located on the ECM fibers for migration. It is assumed that the diffusion of proliferating and quiescent cell population are equal. Therefore, we use a non-linear diffusion function $\Phi_q(\kappa, p, q, v)\nabla q$ for quiescent cells, similar to the diffusion term describing proliferating cells ^[14,15]. The degradation of ECM fibers by proteolytic enzymes creates spatial gradients which directs the migrating cancer cells. Therefore, the function $\Psi(\kappa, v)p\nabla v$ describes the haptotactic movement of cancer cells when they bind to the adhesion sites on the ECM ^[16]. Proliferating cancer cells compete for locally available space with the ECM fibers. This is described by the last but one term in the first equation (2.3) with a proliferating rate $g(c)$ that depends on the concentration $c(x, t)$ of available nutrients. When proliferating cells become deficient in nutrients, they move into a quiescent state. This is modeled by the term $f(c)p$. $f(c)$ describes the rate at which proliferating cells transition into the quiescent state. The function $h(c)$ describes the rate at which quiescent cells die due to nutrients starvation. See the explicit form of the functions $f(c)$ and $h(c)$ as given in ^[15], below. In the model, ECM is non-diffusive, we assume that proliferating and quiescent cells come into contact with the ECM to degrade it at a rate $\partial_v$. The ECM fibers reestablish themselves at a rate $\mu_v$, whiles competing for locally available space with the cells. The competition is modeled by a term similar to that of the proliferating cells. Necrotic cells do not diffuse, rather they are produced when quiescent cells die. A contractivity function $\kappa$ is used to couple the subcellular, microscopic scale biochemical events to the macroscopic scale (population of cells) events. This function determines the strength of tumor cell motility. The contractivity function has a constant time lag $(\tau)$ to account for the delay in the transfer of signals from the subcellular events to affect the dynamics of motile tumor cells on the macroscale. The complete system of equations describing the interactions between proliferating cells, quiescent cells, necrotic cells, ECM, proteolytic enzymes, contractivity function and cell nutrients are given by:

where $\partial_t c = \dfrac{\partial c}{\partial t}$ for the PDEs, $\Omega \in R^d$ for $d = 1$ representing a bounded domain with a smooth boundary. We supplement the model in (2.3) with boundary and initial conditions, in sections 2.3 and 2.4. Whereas $g(c)$ is an increasing function, $f(c)$ and $h(c)$ are both decreasing functions. The rate at which proliferating cells become quiescent is greater than the rate at which quiescent cells die to become necrotic and therefore, $f(c) > h(c)$. $K_r$ and $K_v$ are carrying capacities of cancer cells and ECM respectively, $\delta_u$ is the rate at which cells become necrotic, $\alpha$, $k_1$ and $k_2$ are constants, and $\tau > 0$ is a fixed constant time delay. Additionally, $g(c) = \beta\exp^{\beta c}$, $f(c) = 0.5(1-\tanh(4c-2))$ and $h(c) = \frac{f(c)}{2}$. The non-linear density-dependent diffusion and haptotaxis functions are given by:

$\Phi_p(\kappa, p, q, v) = \frac{D_p\kappa}{1+(\frac{p}{p+q})\frac{v}{K_v}}$, $\Phi_q(\kappa, p, q, v) = \frac{D_q\kappa}{1+(\frac{q}{p+q})\frac{v}{K_v}}$,

$\Psi(\kappa, v) = \frac{D_h\kappa v}{K_v+v}$, $r_3 = (\frac{p}{K_r}+\frac{q}{K_r} +\frac{u}{K_r})$,

2.3.   Boundary conditions

where $\upsilon$ denote the outward unit normal vector on $\partial \Omega$.

2.4.   Initial conditions

$p(0, x) = p_0(x)$, $q(0, x) = q_0(x)$, $u(0, x) = 0$, $\kappa(0, x) = \kappa_0(x)$,

$c(0, x) = c_0(x)$, $y(0, x) = y_0(x)$, $\qquad \qquad \qquad\quad\qquad \qquad \qquad$ $t\in(-\infty, 0]$, $x\in \Omega$

3.   Nondimensionalization

We rescale the system of equations in (2.3) using $p^*: = \frac{p}{K_r}$, $q^*: = \frac{q}{K_r}$, $v^*: = \frac{v}{K_v}$, $u^*: = \frac{u}{K_r}$, $c^*: = \frac{c}{C_0}$, $x^*: = \frac{x}{L}$, $t^*: = \frac{t}{T}$, $t^*: = \eta \theta$ and $y^*: = \frac{y}{l_0}$. $L$ and $T$ are reference length and time scales respectively. We obtain the nondimensionalized system of Eq. (3.1), where the asterix has been dropped for convenience.

where

$D^{*}_p = \frac{D_pT}{L^2}$, $\quad$ $D^{*}_{h} = \frac{D_hTK_v}{L^2}$, $\quad$ $D^{*}_q = \frac{D_qT}{L^2}$, $\quad$ $\lambda = \frac{1}{T}$.

$\delta^{*}_v = \delta_vTK_r$, $\quad$ $\mu^{*}_v = \mu_vT$, $\quad$

$\delta^{*}_u = \delta_u T$

$D^{*}_c = \frac{D_cT}{C_0 L^2}$, $\quad$ $k^{*}_1 = k_1T$, $\quad$ $k^{*}_2 = k_2K_rT$

$k^{+*} = k^{+}\eta TK_v$, $\quad$ $k^{-*} = k^{-}\eta T$

$a^{*} = a\eta T$, $\quad$ $M^{*} = M\eta T$, $\quad$ $\tau^{*} = \frac{\tau}{\eta T}$.

4.   Discretization and numerical experiment

The equations are discretized using the implicit-explicit finite difference method. The space interval $[0, 1]$ is divided into $k$ parts giving us $k+1$ nodes. The space step $\Delta x = 0.01$. We use forward difference for the time derivatives and central difference for the diffusion and haptotaxis terms. Eg. $\dfrac{\partial u}{\partial t} = \dfrac{u^{n+1}-u^{n}}{\Delta t}$, $\dfrac{\partial u}{\partial x} = \dfrac{u_{i+1}-u_{i-1}}{2\Delta x}$ and $\dfrac{\partial^2 u}{\partial x^2} = \dfrac{u_{i+1}-2u_{i}+u_{i-1}}{(\Delta x)^2}$. To simplify the model for discretization purposes we define the following variables: $S_p = \frac{p}{p+q}$, $S_q = \frac{q}{p+q}$, $M = \frac{\kappa}{1+S_jv}$ and $N = \frac{\kappa v}{1+v}p$, where $j\in \{p, q\}$. The evolution equation describing the dynamics of the proliferating cell from equation (3.1) then reads:

Systematically, the system of equations are discretized and numerically solved. It is important to note that the micro-level events update the contractivity function. As mentioned earlier, there is no spatial movement in the equation for ECM density therefore, the ODE is discretized using forward difference in time as follows:

where $n$ denote the time level, $\Delta t = 0.01$ is the time step for the microscale events and $Dt = \eta\Delta t$ is the time step for the macroscale events with $\eta$ as the time scaling factor. The equation describing proteolytic enzyme interactions with ECM and the contractivity functions are also discretized respectively as:

where $k^{+} = 2$, $k^{-} = 0.6$ and $a = 3$. In discretizing the PDE for cancer cells, we use the nonstandard finite difference scheme employed in ^[3]. By this approach, the diffusion and haptotaxis coefficients of the equation for the proliferating cells are handled explicitly while the remaining terms are handled implicitly. The discretization is given below

where ${s_p}^{n+1}_{i+1} = \frac{p^{n+1}_{i+1}}{p^{n+1}_{i+1}+q^{n+1}_{i+1}}$, $\quad$ ${s_p}^{n+1}_{i-1} = \frac{p^{n+1}_{i-1}}{p^{n+1}_{i-1}+q^{n+1}_{i-1}}$, $\quad$ ${s_p}^{n}_{i} = \frac{p^{n}_{i}}{p^{n}_{i}+q^{n}_{i}}$,

$M^{n+1}_{i+1} = \dfrac{\kappa^{n+1}_{i+1}}{1+{s_p}^{n+1}_{i+1}v^{n+1}_{i+1}}$, $\quad$ $M^{n+1}_{i-1} = \dfrac{\kappa^{n+1}_{i-1}}{1+{s_p}^{n+1}_{i-1}v^{n+1}_{i-1}}$, $\quad$ $M^{n}_{i} = \dfrac{\kappa^{n}_{i}}{1+{s_p}^{n}_{i}v^{n}_{i}}$,

$N^{n+1}_{i+1} = \dfrac{\kappa^{n+1}_{i+1}v^{n+1}_{i+1}}{1+v^{n+1}_{i+1}}p^{n+1}_{i+1}$, $\quad$ $N^{n+1}_{i-1} = \dfrac{\kappa^{n+1}_{i-1}v^{n+1}_{i-1}}{1+v^{n+1}_{i-1}}p^{n+1}_{i-1}$ $\quad$, $N^{n}_{i} = \dfrac{\kappa^{n}_{i}v^{n}_{i}}{1+v^{n}_{i}}p^{n}_{i}$.

The index set $R_i = \{i-1, i+1\}$ points at the direct neighbors of the space node $x_i$. The discretization of the equation describing the proliferating cancer cells can thus be written as:

The tridiagonal matrix $\mathbf{A_p}$ has size $(k+ 1)\times (k+1)$, $\mathbf{p^{n+1}}$ is the vector containing the values of $\mathbf{p}$ for the space node $k+1$ at the $(n+1)th$ time level. The vector $\mathbf{\vartheta^n_p}$ contains the entries $Dt \lambda g(c^n_i)p^{n}_i(1-r^n_i-v^{n+1}_i)$ for $i = 1, 2, 3, ..., k$. Similarly, the equation for the quiescent cells read:

where $\chi^{n+1}_{i+1} = \dfrac{\kappa^{n+1}_{i+1}}{1+{s_q}^{n+1}_{i+1}v^{n+1}_{i+1}}$, $\quad$ $\chi^{n+1}_{i-1} = \dfrac{\kappa^{n+1}_{i-1}}{1+{s_q}^{n+1}_{i-1}v^{n+1}_{i-1}}$, $\quad$ $\chi^{n}_{i} = \dfrac{\kappa^{n}_{i}}{1+{s_q}^{n}_{i}v^{n}_{i}}$.

Equation (4.3) can also be written as

The tridiagonal matrix $\mathbf{A_q}$ has size $(k+ 1)\times (k+1)$, $\mathbf{q^{n+1}}$ is the vector containing the values of $\mathbf{q}$ for the space node $k+1$ at the $(n+1)th$ time level. The vector $\mathbf{\vartheta^n_q}$ contains the entries $Dtf(c^n_i)p^n_i$ for $i = 1, 2, 3, ..., k$. From the necrotic cell density and the nutrient concentration we obtain:

where $\mathbf{\vartheta^n_c}$ is the vector containing $\Delta t k_2c_0(1-\alpha(p^n_i+q^n_i))$ for $i = 1, 2, 3, ..., k$.

The following initial conditions were used in the numerical simulations for proliferating cell, quiescent cell, ECM, necrotic cell, a proteolytic enzyme, and contractivity densities respectively.

As mentioned above, early detection of tumor cells is difficult therefore, we initially assume a small lump of tumor cells have penetrated the tissue.

Table 1 shows the parameter description, range, and values used in the model.

Table 1.  Table of parameters used in the model.

Parameter Description	Value	Range	Source
$D_c$ (Diffsion coefficient of nutrient)	$1$	$10^{-2}-2$	[15]
$D_p$ (Diffusion coefficient of proliferating cells)	$10^{-2}$	$10^{-4}-10^{-2}$	[14]
$D_q$ (Diffusion coefficient of quiescent cells)	$10^{-2}$	$10^{-4}-10^{-2}$	[14]
$D_h$ (Haptotaxis coefficient)	$1$	$10^{-3}-10$	[6,14]
$\delta_v$ (Rate of ECM degradation)	$15$	$1-20$	[6]
$a$ (decay rate)	$3$	$1-5$	[6]
$\mu_v$ (Rate of ECM proliferation)	$0.3$	$0.15-2.5$	[3,6]
$c_0$ (Initial concentration of nutrients)	$1$	$0-1$	[15]

---

 | Show Table

DownLoad: CSV

5.   Discussion

In this section, we show the numerical solutions for the densities of three cell layers (proliferating, quiescent, necrotic), ECM, and contractivity function in space and specific intervals of time ($t = 1, 2, 5, 10, ...50$). Figure 1 (A-L), Figure 2 (A-F), and Figure 3 (A-C) show qualitative solutions without delay ($\kappa = 1$, solid lines) and with a delay ($\tau = 5$ dashed lines) for: $\eta = 0.01, \lambda = 2$; $\eta = 0.01, \lambda = 0.02$: and $\eta = 0.005, \lambda = 1.5$ respectively.

Figure 1.  Evolution of proliferating tumor cell density (blue), quiescent cell density (pink), necrotic cell density (red), ECM fiber density (green), and contractivity function (black) in the pure macroscopic setting (solid lines) and with the delay of $\tau = 5$ (dashed line), $\eta = 0.01$, $M = 1$, $\lambda = 2$ for equations (3.1) and (2.4).

DownLoad: Full-Size Img PowerPoint

Figure 2.  Evolution of proliferating tumor cell density (blue), quiescent cell density (pink), necrotic cell density (red), ECM fiber density (green), and contractivity function (black) in the pure macroscopic setting (solid lines) and with the delay of $\tau = 5$ (dashed line), $\eta = 0.01$, $M = 1$, $\lambda = 0.02$ for equations (3.1) and (2.4).

DownLoad: Full-Size Img PowerPoint

Figure 3.  Evolution of proliferating tumor cell density (blue), quiescent cell density (pink), necrotic cell density (red), ECM fiber density (green), and contractivity function (black) in the pure macroscopic setting (solid lines) and with the delay of $\tau = 5$ (dashed line), $\eta = 0.005$, $M = 1$, $\lambda = 1.5$ for equations (3.1) and (2.4).

DownLoad: Full-Size Img PowerPoint

As earlier mentioned, biochemical reaction at the subcellular level takes place on a short time scale whiles its effects on the motile behavior of tumor cells take a much longer time to happen. We coupled the two-time scales by introducing the delay term in the contractivity function which we modeled as an input to the nonlinear density-dependent diffusion function of proliferating and quiescent cells, as well as the haptotaxis coefficient of proliferating cells. As time progress, proliferating cells invade faster without the delay ($\kappa = 1$) than with delay as expected. Therefore it can be concluded that subcellular events slow down the invasion of proliferating cells into the tissue. Furthermore, there are higher peaks at the leading edge for proliferating cells (this may be due to the cells encountering undegraded tissue or bone).

The density of necrotic cells builds up gradually at the center ($x = 0$) as time increases. Accounting for the subcellular events shows that they seem to have no effect on the evolution of quiescent and necrotic cell densities. The invasion of ECM fiber is as expected. From Figures 1, 2 and 3, we notice that the solution has not reached the final time. Therefore, there is negligible influence of the boundary conditions on the invasion patterns.

6.   Conclusion

The study makes the following contributions to existing mathematical models on tumor invasion of host tissues: Our numerical solution for the system of equations showed the three distinct layers of proliferating, quiescent and necrotic cell densities. The subcellular, microscopic scale and macroscale events have varying spatial and time scales, we coupled the two-time scales using a contractivity function. As a result, we were able to determine the effects of biochemical dynamics at the subcellular level on the strength and motility of tumor cells. The model accounts for contact inhibition of proliferating and quiescent cells with the tumor micro-environments.

Conflict of interest

The authors declare there is no conflicts of interest in this paper.