Mathematical Modeling of Dynamics of Cancer Invasion in Human Body Tissues

Preprint

Review

Mathematical Modeling of Dynamics of Cancer Invasion in Human Body Tissues

Altmetrics

Downloads

211

Views

Comments

Supplementary Material

supplementary.pdf (695.31KB )

Dinnah Mogire^*,

Joash Kerongo,Vincent Bulinda

Dinnah Mogire^*,

Joash Kerongo,Vincent Bulinda

This version is not peer-reviewed

This preprints belongs to the Topic

Mathematical Modeling

Submitted:

09 April 2023

Posted:

10 April 2023

Read the latest preprint version here

Alerts

Abstract

Cancer diseases lead to the second-highest death rate all over the world. The dynamics of invasion of cancer cells into the human body tissues and metastasis are the main causes of death in patients with cancer. This study deals with theoretical investigation of the dynamics of invasion of cancer cells for tumour growths in human body tissues using discretized Cahn-Hilliard, concentration and reaction-diffusion equations which were solved by Finite Difference Method with the aid of MATLAB computer software. A Crank-Nicolson numerical scheme was developed for the discretized model equations. The numerical result obtained was used to describe the dynamics of cancer invasion of tissues with respect to cancer cells density on tumour growth, turbulence and mobility and equilibrium between charge and discharge of cancer cells. The results of the study provide new insights into combating cancer disease by providing mitigating and intervention measures to this major health problem.

Keywords:

Subject: Computer Science and Mathematics - Applied Mathematics

Nomanclature:

Μ

, Mobility Parameter (m²/V.\ s);

N

, Number of living cells per unit volume in Nanometers (

n m

); S_T, Net sources of tumour cells (nm); S_D_, Net sources of dead cells (nm);

μ

, Density of the tumour (kg (m^-3);

F (ϕ)

Double well bulk energy(joule);

ε

, Parameter related to the thickness of thediffuse interface tumour and host domains;

u

, Tumour velocity(m/s);

β (t) \sqrt{F (ϕ)}

; Space–time dependent Lagrange multiplier;

p

, Pressure (N/m²);

n

, Nutrient concentration;

ϕ

, Sum of the volume fractions of viable and dead tumour cells;

ε

, Host tissue;

ψ

, Dead tumour cells;

P

, Population density;

v

, Advection cancer invasion velocity;

D

, Cancer cells diffusivity (

μ M

²/h);

ρ

, Density of a single tumour cell (

k g . m^{3}

);

D_{m}

, MDE diffusion coefficient;

β

, MDE production rate by cells ;

γ

, MDE natural decay rate;

Ω (t)

, Area occupied by tumour at time

t

;

C_{p}

, Specific heat capacity (

J . k g^{- 1} . k^{- 1}

);

α,

Repose angle;

κ

, Materials conductivity

(W . m^{- 1} . k^{- 1})

;

δ

, Flowing layer depth;

\Pr

, Prandtl number;

P e

, Peclet number;

Re

, Reynolds number;

S c

, Schmidt number;

J_{x}

Total flux of the cancer cells

x

-axis;

J_{y}

, Total flux of the cancer cells along the

y

-axis.

1. Introduction

Cancer remains the world's primary cause of death, despite continuous therapeutic advances, making it a major public health issue with devastating societal consequences (Perthame, 2016 and Zhang et al., 2022). From what can be gleaned from the databases of the World Health Organization (WHO) and the United Nations (UN), it is clear that cancer is currently, and is likely to remain, one of the top causes of mortality around the world. Cancer ranks high on the list of global killers and causes a great deal of suffering. Annually, it's responsible for the deaths of around six million individuals. In light of this, it is clear why efforts toward cancer prevention, treatment, and discovery are so vital (Kolev and Zubik-Kowal, 2010). The anticipated 18.1 million new cases and 9.6 million deaths from cancer worldwide in 2018 demonstrate the continued escalation in this disease's global burden. It's putting a major burden on healthcare systems and populations worldwide. WHO's Global Cancer Observatory projects that ten million people will lose their lives to cancer in 2020, up from nine million in 2018, and that number will climb to over sixteen and a half million by 2040 if preventative measures aren't implemented (WHO, 2020 & 2018, Roose et al, 2007).

Some of the biggest obstacles to a comprehensive understanding of cancer development and treatment are metastatic colonization, dormancy, relapse, multiple drug therapy, immune resistance, and the use of mathematical models and simulations to predict the course of the disease and the optimal treatment. What makes this phenomenon so difficult to understand is that cancer is an ancient disease that appears to be inherent to the complexity of creatures produced by evolution, in contrast to viral infections, which only became major killers with the rise of vast societies a few hundred years ago (Perthame, 2016).

Cancer is a disease characterized by a series of evolutionary steps at the somatic level, including tumour initiation, development, dissemination, immune evasion, and the emergence and maintenance of drug resistance. The evolutionary history of a tumour can be deduced from molecular data, and evolutionary theory can be used to examine the dynamics of tumour cell populations. The promise of several approaches to modeling cancer's development has been recently assessed. They include the phylogenetic approaches to modeling the evolutionary relationship between tumor subclones and the probabilistic graphical models for characterizing dependencies among mutations, as well as models of tumour onset and progression based on the dynamics of the population. Understanding the genesis of tumours, as well as the likely course of a disease and the effectiveness of treatments like targeted therapy, can be greatly aided by evolutionary modeling. Hence, the multistep theory is incorporated into the population dynamics models. Tumour phylogenetic tree, single-cell approach, hybrid models, and stochastic models are some of the phylogenetic approaches (Beerenwinkel, Schwarz, Gerstung and Markowetz 2015).

In order to deal with the difficulties associated with tumour growth, Perthame (2016) provides some invasion models. First, there are models of interplanetary invasion based on the Reaction-Diffusion Equations. Considered alongside constrained expansion and dispersal via diffusion, it depicts the spread of an infectious disease, the arrival of a new species (with no natural enemies), a combustion wave, or a new way of thinking. Second, the Fisher-KPP model is suggested to model the spatiotemporal development of a tumour under chemotherapeutic influence.

In order to tackle the many difficulties associated with tumour development, several models based on Ordinary Differential Equations (ODEs) have been presented. Two species, one predator and one prey, interacting in a biological system is the focus of several theoretical frameworks, such as the Lotka-Volterra model. cells that are both active and dormant, or P-Q, Model; According to the NCI's Dictionary of Cancer Terminology, cell proliferation refers to the increase in cell population that occurs as a result of normal cellular processes including growth and division. An absence of proliferation that is both transitory and reversible characterizes a state called quiescence. To wit. When it comes to describing and analyzing the immune response to cancer growth, the ODEs models are invaluable tools (Fornier & Sagot, 2011).

Tumour growth models have been demonstrated to be an essential tool in the development of an engineering foundation for cancer therapy through the use of therapeutic procedure design coupled with control engineering or the use of the models for simulation and assessment of treatment processes. This work further demonstrates the usefulness of mathematical modeling in describing tumour formation, a process of great complexity, by highlighting the characteristics and physiological processes of the tumours. For this reason, they built a partial differential equation-based mathematical model for tumour growth and presented its results in their study. In this model, the densities of proliferative, quiescent, necrotic, and surrounding cells are described together with the flow of nutrients using a deterministic model of an avascular tumour growth framed in a system of nonlinear coupled PDEs (Akhtar et al., 2021).

This paper is organized as follows. First, the background to the study is described in section1. Then, the governing mathematical principles and equations are given in section

2

, model formulation and applications are presented in section

3

, results and discussion on results is provided in section 4, summary and conclusion in section

5

, and recommendations in section

6

2. Background

The mechanisms involved in the evolution of cancer are extremely intricate and dynamic. Due to its enormous complexity, cancer has become a very effective human killer in modern times by hindering early identification and treatment. The numerous methods and approaches used to treat cancer aim to stop the disease's spread at different phases (such as immortalization or unending cell development, transformation, and metastasis). Yet, the bigger achievements in science and medicine have not yet connected all the dots, leaving us without an apparent solution at this time. We are still unable to develop a vaccine or an effective anti-cancer medicine since we still don't fully understand how cancer develops and manifests (Fornier & Sagot, 2011).

One such area that pertains to the cutting-edge study of cancer growth and treatment is mathematical modeling. The goal of developing mathematical models of cancer is to forecast tumour growth and treatment options. A highly associative topic, mathematical modeling integrates the real-world issue with multiple simulative models built using mathematical techniques.

There is still tremendous space for improvement in cancer therapy, despite the fact that it has been more effective for many tumours in recent years. Knowing the basics of how tumour cells invade can help doctors prescribe treatments and preventative measures for cancer. This cancer invasion process can be mathematically described as a system (or systems) of differential equations, which, when solved, will essentially offer a useful contribution to solving the cancer problem. In order to maximize the cancer problem, mathematical oncology will be employed to deliver an integrated solution without the need for clinical trials and experimentation. Therefore, there will be need to determine integrable theoretical solution to the dynamics of cancer invasion in reference to density of tumour cells, turbulence and motility of tumour cells and to determine the equilibrium between charge and discharge of tumour cells. (Socolofsky, S. A. and Jirka.G.H., 2002).

2.1. Governing Equations and Theories

This section discusses the equations that govern the growth rate, concentration and invasion of cancer cells with reference to its density, growth rate, degradation, renewal, mobility and concentration in the tumour. Since this is a two-dimensional problem, the equations are presented along the horizontal and vertical axes. The following universal laws form the basis for the general equations in the dynamics of cancer cells: The equations of diffusion and concentration, as well as the Cahn-Hilliard equation. In order to facilitate their incorporation into the aforementioned computer code, these essential equations are presented in Cartesian form.

It is not possible to get analytical solutions to the equations driving cancer invasion and tumour growth since they are extremely non-linear PDEs. The Navier-Stokes equations are used to characterize the dynamics of fluids. Hence, the Navier-Stokes equations are used to model the flow of fluids, such as the weather, ocean currents, water flow in a pipe, and air flow over a wing. Its use as models extends to other fields, such as the analysis of pollutants, the development of power plants, and the research of blood flow. These equations include the following;

2.1.1. Continuity Equation

\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0

(2.1)

Along with a mathematical definition of the relevant process, the equation for conservation of mass (or continuity equation) can be used to derive a differential equation describing the movement of cancer cells (Sarah, 2011 & Elaine, 2009). The conservation equation (2.1) can be examined under two distinct circumstances. The first scenario is satisfied when the charging and discharging rates of the cancer cells are equal. The second scenario involves investigating what occurs when the charging and discharging rates of cancer cells are different. The first condition will be considered in this study to develop the dynamics of cancer cells in human body.

2.1.2. Momentum Equation

According to the principle of momentum conservation in fluid dynamics, the momentum of a closed system does not change. The momentum equation, a form of Newton's Second Law, describes the relationship between the total force acting on an element and the rate of change in that element's momentum or acceleration. These partial differential equations (PDEs) provide the incompressible version of the Navier-Stokes equations. The Navier-Stokes equation in x-direction is given as

ρ (\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + u \frac{\partial u}{\partial y}) = - \frac{\partial p}{\partial x} + μ (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}}) + ρ f

(2.2)

The Navier-Stokes equation in y-direction is given as

ρ (\frac{\partial v}{\partial t} + v \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}) = - \frac{\partial p}{\partial y} + μ (\frac{\partial^{2} v}{\partial x^{2}} + \frac{\partial^{2} v}{\partial y^{2}}) + ρ f

(2.3)

2.1.3. Concentration Equation

In fluid mechanics, the concentration equation (derived from the scalar transport equation) describes how material is transported (and, thus, mixed) from a continuum viewpoint. In this equation, material (with concentration

C)

is passively carried by a velocity field

u

while simultaneously undergoing diffusion (where

D

is the diffusion coefficient). In contrast to the Navier-Stokes equations, the velocity field

u

is known, and with the given initial condition, the evolution of the concentration can be determined for all time. The application of the concentration equation to physical systems can shed light onto how these systems diffuse for different parameters.

The net movement of the cancer cells in

x

and

y

respectively are

\begin{array}{l} u δ t \pm δ x \\ u δ t \pm δ y \end{array}

(2.4)

Thus, the total flux in the

x

and

y

directions are

J_{x}

and

J_{y}

, respectively including the advective transport and a Fickian diffusion term, we obtain

J_{x} = u C + C_{x} = u C - D \frac{\partial C}{\partial x}

(2.5)

J_{y} = v C + C_{y} = v C - D \frac{\partial C}{\partial y}

(2.6)

\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = D (\frac{\partial^{2} C}{\partial x^{2}} + \frac{\partial^{2} C}{\partial y^{2}})

(2.7)

where

u (x, y, t)

is the given divergence-free advective velocity vector field,

t

is time and

C_{(x, y, t)}

is the fluid concentration levels in mixing process,

C (x, y, t)

is the concentration of cancer cells and

D

is the diffusion coefficient.

2.1.4. Energy Equation

Energy can neither be created nor be destroyed; can only change physical form. Energy equation can be written in many different ways, such as the one given below

ρ C_{P} (\frac{\partial θ}{\partial t} + u \frac{\partial θ}{\partial x} + v \frac{\partial θ}{\partial y}) = κ \nabla^{2} θ + μ ϕ

(2.8)

2.1.5. Cahn–Hilliard Equation

The original governing equations for the tumour growth in Wise et al., (2008), are the Cahn–Hilliard equation

\frac{\partial ϕ}{\partial t} = M ρ (\frac{\partial^{2} ϕ}{\partial x^{2}} + \frac{\partial^{2} ϕ}{\partial y^{2}} + \frac{\partial^{2} ϕ}{\partial z^{2}}) - v (\frac{\partial ϕ}{\partial x} + \frac{\partial ϕ}{\partial y} + \frac{\partial ϕ}{\partial z}) + S_{T} ϕ

(2.9)

μ = (F^{'} (ϕ) - ε^{2} Δ ϕ)

(2.10)

u = \nabla P - \frac{γ}{ε} ϕ \nabla μ

(2.11)

where, mobility

M > 0

The diffuse interface model of Wise et al., (2008) is reformulated using the conservative second-order AC equation with a space–time dependent Lagrange multiplier in Equation (2.11).

The diffuse interface model of Wise et al (2008), involves the second-order CH equation with a source term. The proposed model consists of the following three equations

\frac{\partial ϕ}{\partial t} = M ρ (\frac{\partial^{2} ϕ}{\partial x^{2}} + \frac{\partial^{2} ϕ}{\partial y^{2}} + \frac{\partial^{2} ϕ}{\partial z^{2}}) - v (\frac{\partial ϕ}{\partial x} + \frac{\partial ϕ}{\partial y} + \frac{\partial ϕ}{\partial z}) + S_{T} ϕ

(2.12)

\frac{\partial ψ}{\partial t} = M μ (\frac{\partial^{2} ψ}{\partial x^{2}} + \frac{\partial^{2} ψ}{\partial y^{2}}) - u (\frac{\partial ψ}{\partial x} + \frac{\partial ψ}{\partial y}) + ϕ S_{D}

(2.13)

\frac{\partial ϕ}{\partial t} = M ϕ (- F^{'} (ϕ) + c^{2} (\frac{\partial ϕ}{\partial x} + \frac{\partial ϕ}{\partial y}) β (t) F (ϕ))

(2.14)

3. Mathematical Formulation

Our model will be solved numerically in two dimensions. With the help of the finite difference approximations, we conduct a spatial discretization of the model equations, reducing them to a set of (time-dependent) partial differential equations that can be solved with the help of the central difference formula.

In Figure 3.1, x is the horizontal coordinate and y is the vertical coordinate. The flowing layer depth is

δ

, and the length of the tissue under consideration is

L

Figure 3.1. Schematic cancer cell flow chamber.

There is a cylinder-shaped flow chamber. Cancer cells, which are seen as small, circular particles, and blood, which is depicted as the fluid represented by the arrow, enter the chamber and become well mixed. A current of fluid is moving through the chamber. Arteries, veins, and capillaries are the blood vessels (tubes) via which blood travels throughout the body. The cardiovascular system functions like this. It is used to transport gases such as oxygen, carbon dioxide, and others. The presence of blood supply is crucial in the metastatic process. Some cancer cells are able to leave the main tumour and enter the bloodstream by squeezing through tiny crevices in the walls of blood vessels. Modern cancer cells in the circulation are able to adhere to the walls of blood vessels, creating microscopic channels through which they can inject genetic information that changes the endothelial cells lining the blood vessels, making them far more friendly to other cancer cells.

In nature, transport occurs in fluids through the combination of advection and diffusion. This section incorporates advection into our concentration equation.

Figure 3.2. Schematic diagram of a control volume with cross flow.

The concentration equation is derived using the principle of superposition, which states that linearly independent processes such as advection and diffusion can be put together. Each cancer cell, over the course of time, will diffuse to the left or right by a single step. Due to diffusion, each cancer cells in time

t

will move either one step to the left or one step to the right ( ±

δ x

). Due to advection, each cancer cells will also move

u δ t

in the cross-flow direction. The crossflow does not affect whether or not cancer cells diffuse to the right or left. However, its presence only amplifies the effect of the original process.

3.4. Methods of Solution

The Finite Difference Technique was used for this research. Cahn-Hilliard, diffusion, and concentration equations are solved by means of two numerical schemes: a central Crank-Nicolson scheme and a central Difference scheme. By discretizing the supplied equation and developing numerical techniques analogous to the equations, the approaches obtain a finite system of linear or nonlinear algebraic equations from the Partial Differential Equations. Using the specified boundary conditions in mind, we solve the equations using MATLAB software.

3.4.1. Central Difference Scheme

The Central Difference scheme is a finite difference method used to get numerical solutions to a differential equation in applied mathematics by optimizing the approximation for the differential operator at the center of the patch under discussion. Its benefits include a quicker convergence rate than other finite differencing methods like the forward and backward differencing methods, as well as being straightforward to comprehend and apply at least for simple material relations.

3.4.2. Crank Nicolson Scheme

In a central Crank Nicolson scheme, the temporal variable

t

is replaced by forward difference scheme:

y

is replaced by central difference scheme while

x

is by the average of central difference approximation at

j

th and (

j + 1

)th level.

3.5. Dimensionalizing Tumour Concentration Equation

Our goal is to solve the concentration equation for tumour concentration

C (x, y, t)

. In dimensionless variables, the governing equation, Landau and Lifshitz (1959) is;

\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = \frac{1}{P e} (\frac{\partial^{2} C}{\partial x^{2}} + \frac{\partial^{2} C}{\partial y^{2}}),

(3.1)

where

u (x, y, t)

and

v (x, y, t)

are the given divergence-free advective velocity vector fields along

x

and

y

directions respectively ,

t

is time and

C (x, y, t)

is the tumour concentration levels in human body tissues.

But the P’eclet number is given by

P e = \frac{U L}{D}

(3.2)

Pe is the P´eclet number, where

U

is an axial velocity,

L

is arterial length, and

D

is the diffusivity of the tumour cells.

Substituting (3.2) into (3.1) gives

\frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = P e (\frac{\partial^{2} C}{\partial x^{2}} + \frac{\partial^{2} C}{\partial y^{2}})

(3.3)

3.6. Discretization of Concentration Equation

We discretize the concentration partial differential equation and form a central Crank Nicolson scheme which we eventually solve using the finite difference method. In application, equation (3.3) is discretized to study the effects of Peclet number for tumour concentration levels. Using a central difference numerical scheme,

C_{t}

is replaced by forward difference scheme while

C_{x x}

is the average of central difference approximation at

j

th and

j + 1

th level, and

C_{y y}

is replaced by central difference approximation .When these approximations are substituted into equation (3.3), and taking

u = v = 1

, the second derivative in concentration equation at node could be represented as follows:

\frac{C_{i, j}^{n + 1} - C^{n}_{i, j}}{Δ t} + \frac{C_{i + 1, j}^{n} - C^{n}_{i - 1, j}}{2 Δ x} + \frac{C_{i, j + 1}^{n} - C^{n}_{i, j - 1}}{2 Δ y} = P e [\frac{1}{2} (\frac{C_{i + 1, j}^{n} - 2 C_{i, j}^{n} + C_{i - 1, j}^{n}}{{(Δ x)}^{2}} + \frac{C_{i + 1, j + 1}^{n} - 2 C_{i, j + 1}^{n} + C_{i - 1, j + 1}^{n}}{{(Δ x)}^{2}}) + \frac{C_{i, j + 1}^{n} - 2 C_{i, j}^{n} + C_{i, j - 1}^{n}}{{(Δ y)}^{2}}]

(3.4)

We investigate the effect of Pe number on the tumour concentration levels. Taking

ϕ = \frac{Δ t}{(Δ x)} = \frac{Δ t}{(Δ y)}

, and

μ = \frac{Δ t}{{(Δ x)}^{2}} = \frac{Δ t}{{(Δ y)}^{2}}

Δ x = Δ y

and multiply by

2 Δ t

throughout (3.4) and re-arranging, we get the scheme

(ϕ - μ P e) C_{i + 1, j}^{n} - (ϕ + μ P e) C_{i - 1, j}^{n} + (6 μ P e - 2) C_{i, j}^{n} = - ϕ C_{i, j + 1}^{n} + μ P e C_{i - 1, j + 1}^{n} + (2 μ P e + ϕ) C_{i, j - 1}^{n} + μ P e C_{i + 1, j + 1}^{n} - 2 C_{i, j}^{n + 1}

(3.5)

Taking

Δ x = Δ y = 0.1

, and

Δ t = 0.01

\Rightarrow ϕ = 0.1

and

μ = 2

we obtain the central difference scheme as

(0.1 - 2 P e) C_{i + 1, j}^{n} - (0.1 + 2 P e) C_{i - 1, j}^{n} + (12 P e - 2) C_{i, j}^{n} = - 0.1 C_{i, j + 1}^{n} + 2 P e C_{i - 1, j + 1}^{n} + (4 P e + 0.1) C_{i, j - 1}^{n} + 2 P e C_{i + 1, j + 1}^{n} - 2 C_{i, j}^{n + 1}

(3.6)

The above central difference scheme can be in form of six algebraic equations when i = 1, 2, 3,4,5……

j = 1 and n = 0 as

\begin{array}{l} (0.1 - 2 P e) C_{2, 1}^{0} - (0.1 + 2 P e) C_{0, 1}^{0} + (12 P e - 2) C_{1, 1}^{0} = - 0.1 C_{1, 2}^{0} + 2 P e C_{0, 2}^{0} + (4 P e + 0.1) C_{1, 0}^{0} + 2 P e C_{2, 2}^{0} - 2 C_{1, 1}^{1} \\ (0.1 - 2 P e) C_{3, 1}^{0} - (0.1 + 2 P e) C_{1, 1}^{0} + (12 P e - 2) C_{2, 1}^{0} = - 0.1 C_{2, 2}^{0} + 2 P e C_{1, 2}^{0} + (4 P e + 0.1) C_{2, 0}^{0} + 2 P e C_{3, 2}^{0} - 2 C_{2, 1}^{1} \\ (0.1 - 2 P e) C_{4, 1}^{0} - (0.1 + 2 P e) C_{2, 1}^{0} + (12 P e - 2) C_{3, 1}^{0} = - 0.1 C_{3, 2}^{0} + 2 P e C_{2, 2}^{0} + (4 P e + 0.1) C_{3, 0}^{0} + 2 P e C_{4, 2}^{0} - 2 C_{3, 1}^{1} \\ (0.1 - 2 P e) C_{5, 1}^{0} - (0.1 + 2 P e) C_{3, 1}^{0} + (12 P e - 2) C_{4, 1}^{0} = - 0.1 C_{4, 2}^{0} + 2 P e C_{3, 2}^{0} + (4 P e + 0.1) C_{4, 0}^{0} + 2 P e C_{5, 2}^{0} - 2 C_{4, 1}^{1} \\ (0.1 - 2 P e) C_{6, 1}^{0} - (0.1 + 2 P e) C_{4, 1}^{0} + (12 P e - 2) C_{5, 1}^{0} = - 0.1 C_{5, 2}^{0} + 2 P e C_{4, 2}^{0} + (4 P e + 0.1) C_{5, 0}^{0} + 2 P e C_{6, 2}^{0} - 2 C_{5, 1}^{1} \end{array}}

(3.7)

We apply the following initial conditions and boundary conditions

C (x, y, 1) = 0, t = 0

(3.8)

C (0, y, t) = 10, C (x, 0, t) = C (x, 2, t) = 0, t > 0, x \neq y

(3.9)

The above algebraic equations in (3.9) can be written in matrix form as

[\begin{matrix} (12 P e - 2) & (0.1 - 2 P e) & 0 & 0 & 0 & 0 \\ - (0.1 + 2 P e) & (12 P e - 2) & (0.1 - 2 P e) & 0 & 0 & 0 \\ 0 & - (0.1 + 2 P e) & (12 P e - 2) & (0.1 - 2 P e) & 0 & 0 \\ 0 & 0 & - (0.1 + 2 P e) & (12 P e - 2) & (0.1 - 2 P e) & 0 \\ 0 & 0 & 0 & - (0.1 + 2 P e) & (12 P e - 2) & (0.1 - 2 P e) \\ 0 & 0 & 0 & 0 & - (0.1 + 2 P e) & (12 P e - 2) \end{matrix}] [\begin{matrix} C_{1, 1}^{0} \\ C_{2, 1}^{0} \\ C_{3, 1}^{0} \\ C_{4, 1}^{0} \\ C_{5, 1}^{0} \end{matrix}] = [\begin{matrix} 1 + 40 P e \\ 20 P e \\ 20 P e \\ 20 P e \\ 20 P e \end{matrix}]

(3.10)

Solving the matrix system (3.10) using MATLAB, we obtain the solutions for varying values of Pe to get results as in Table 4.3.

3.7. Cahn-Hilliard Model

The original governing equations for the tumour growth in Wise et al (2008), are the Cahn–Hilliard equation

\frac{\partial ϕ}{\partial t} = M ρ (\frac{\partial^{2} ϕ}{\partial x^{2}} + \frac{\partial^{2} ϕ}{\partial y^{2}} + \frac{\partial^{2} ϕ}{\partial z^{2}}) - v (\frac{\partial ϕ}{\partial x} + \frac{\partial ϕ}{\partial y} + \frac{\partial ϕ}{\partial z}) + S_{T} ϕ

(3.11)

where

ρ = (F^{'} (ϕ) - ε^{2} Δ ϕ)

v = \nabla P - \frac{γ}{ε} ϕ \nabla ρ

Where

M > 0

is mobility,

ρ

is the density of tumour cells,

u

is the tumour velocity and S_T is the source of tumour cells.A numerical algorithm using Finite difference method is described in this section.The diffuse interface model of Wise et al (2008), involves the second-order CH equation with a source term. We consider variables

M

and

ρ

, and take effect of

u

and

S_{T}

as negligible and constant and vary

M

and

ρ

. The variables are taken as

M = 10, 15, 20

and

μ = 1000, 2000, 3000

3.8. Discretization of Cahn-Hilliard Model

We discretize the concentration partial differential equation and form a central Crank Nicolson scheme which we eventually solve using the finite difference method. Equation (3.16) is discretized to study the effects of mobility

M

, and density

ρ

, of tumour cells. Using a central difference numerical scheme,

ϕ_{t}

is replaced by forward difference scheme while

ϕ_{x x}

is the average of central difference approximation at

j

th and

j + 1

th level, and

ϕ_{y y}

is replaced by central difference approximation .When these approximations are substituted into equation (3.11), and according to Crank and Nicolson (1947) the second derivative in Cahn–Hilliard equation at nodal point. Letting the superscript

n

indicate the time (temporal variable) while the subscript

i

and

j

indicate the spatial variable

x

and

y

Discretization of equation (3.16) gives

\frac{ϕ_{i, j}^{n + 1} - ϕ_{i, j}^{n}}{Δ t} = Μ ρ (\frac{1}{2} (\frac{(ϕ_{i + 1, j}^{n} - 2 ϕ_{i, j}^{n} + ϕ_{i - 1, j}^{n})}{{(Δ x)}^{2}} + \frac{ϕ_{i + 1, j + 1}^{n} - 2 ϕ_{i, j + 1}^{n} + ϕ_{i - 1, j + 1}^{n}}{{(Δ x)}^{2}}) + \frac{ϕ_{i, j + 1}^{n} - 2 ϕ_{i, j}^{n} + ϕ_{i, j - 1}^{n}}{{(Δ y)}^{2}}) - v (\frac{ϕ_{i + 1, j}^{n} - ϕ_{i - 1, j}^{n}}{(2 Δ x)} + \frac{ϕ_{i, j + 1}^{n} - ϕ_{i, j - 1}^{n}}{(2 Δ y)}) + \frac{ϕ_{i + 1, j}^{n} + ϕ_{i, j}^{n}}{2} S_{T}

(3.12)

Similarly, in three dimension the subscript

i

and

k

indicate the spatial variable

x

and

z

\frac{ϕ_{i, k}^{n + 1} - ϕ_{i, k}^{n}}{Δ t} = Μ ρ (\frac{1}{2} (\frac{(ϕ_{i + 1, k}^{n} - 2 ϕ_{i, k}^{n} + ϕ_{i - 1, k}^{n})}{{(Δ x)}^{2}} + \frac{ϕ_{i + 1, k + 1}^{n} - 2 ϕ_{i, k + 1}^{n} + ϕ_{i - 1, k + 1}^{n}}{{(Δ x)}^{2}}) + \frac{ϕ_{i, k + 1}^{n} - 2 ϕ_{i, k}^{n} + ϕ_{i, k - 1}^{n}}{{(Δ y)}^{2}}) - v (\frac{ϕ_{i + 1, k}^{n} - ϕ_{i - 1, k}^{n}}{(2 Δ x)} + \frac{ϕ_{i, k + 1}^{n} - ϕ_{i, k - 1}^{n}}{(2 Δ y)}) + \frac{ϕ_{i + 1, k} + ϕ_{i, k}^{n}}{2} S_{T}

(3.13)

We will continue with variables

x

and

y

. Take constant values of

v

and

S_{T}

so that

v = S_{T} = 1

and multiply by 2

Δ t

with square mesh for variables

Δ x = Δ y

We have

ϕ_{i, j}^{n + 1} - ϕ_{i, j}^{n} = \frac{M ρ Δ t}{{(Δ x)}^{2}} (ϕ_{i + 1, j}^{n} - 2 ϕ_{i, j}^{n} + ϕ_{i - 1, j}^{n} + ϕ_{i + 1, j + 1}^{n} - 2 ϕ_{i, j + 1}^{n} + ϕ_{i - 1, j + 1}^{n} + 2 ϕ_{i, j + 1}^{n} - 4 ϕ_{i, j}^{n} + 2 ϕ_{i, j - 1}^{n}) - \frac{Δ t}{Δ x} (ϕ_{i + 1, j}^{n} - ϕ_{i - 1, j}^{n} + ϕ_{i, j + 1}^{n} - ϕ_{i, j - 1}^{n}) + Δ t (ϕ_{i + 1, j}^{n} + ϕ_{i, j}^{n})

(3.14)

Take the mesh sizes

r = \frac{Δ t}{{(Δ x)}^{2}}, τ = \frac{Δ t}{Δ x}, ξ = Δ t

\begin{array}{l} ϕ_{i, j}^{n + 1} - ϕ_{i, j}^{n} = M ρ r (ϕ_{i + 1, j}^{n} - 2 ϕ_{i, j}^{n} + ϕ_{i - 1, j}^{n} + ϕ_{i + 1, j + 1}^{n} - 2 ϕ_{i, j + 1}^{n} + ϕ_{i - 1, j + 1}^{n} + 2 ϕ_{i, j + 1}^{n} - 4 ϕ_{i, j}^{n} + 2 ϕ_{i - 1, j - 1}^{n}) \\ - τ (ϕ_{i + 1, j}^{n} - ϕ_{i - 1, j}^{n} + ϕ_{i, j + 1}^{n} - ϕ_{i, j - 1}^{n}) + ξ (ϕ_{i + 1, j}^{n} + ϕ_{i, j}^{n}) \end{array}

(3.15)

But if

Δ t = 0.01, Δ x = 0.1

r = \frac{0.01}{{(0.1)}^{2}} = 1, τ = \frac{0.01}{(0.1)} = 0.1, ξ = \frac{0.01}{1} = 0.01

Equation (3.15) becomes

\begin{matrix} ϕ_{i, j}^{n + 1} - ϕ_{i, j}^{n} = M ρ (ϕ_{i + 1, j}^{n} - 2 ϕ_{i, j}^{n} + ϕ_{i - 1, j}^{n} + ϕ_{i + 1, j + 1}^{n} - 2 ϕ_{i, j + 1}^{n} + ϕ_{i - 1, j + 1}^{n} + 2 ϕ_{i, j + 1}^{n} - 4 ϕ_{i, j}^{n} + 2 ϕ_{i, j - 1}^{n}) \\ - 0.1 (ϕ_{i + 1, j}^{n} - ϕ_{i - 1, j}^{n} + ϕ_{i, j + 1}^{n} - ϕ_{i, j - 1}^{n}) + 0.01 (ϕ_{i + 1, j}^{n} + ϕ_{i, j}^{n}) \end{matrix}

(3.16)

We multiply both sides by 100

100 ϕ_{i, j}^{n + 1} - 100 ϕ_{i, j}^{n} = 100 M ρ (ϕ_{i + 1, j}^{n} - 2 ϕ_{i, j}^{n} + ϕ_{i - 1, j}^{n} + ϕ_{i + 1, j + 1}^{n} - 2 ϕ_{i, j + 1}^{n} - ϕ_{i - 1, j + 1}^{n} + 2 ϕ_{i, j + 1}^{n} - 4 ϕ_{i, j}^{n} + 2 ϕ_{i, j - 1}^{n}) + 100 ϕ_{i + 1, j}^{n} - 100 ϕ_{i - 1, j}^{n} + 100 ϕ_{i, j + 1}^{n} - 100 ϕ_{i, j - 1}^{n} + ϕ_{i + 1, j}^{n} + ϕ_{i, j}^{n}

(3.17)

Putting the

j + 1

level on LHS

- 20 M ρ ϕ^{n}_{i + 1, j + 1} + (1 - 50 M ρ) ϕ^{n}_{i, j + 1} - 20 M ρ ϕ^{n}_{i - 1, j + 1} = (10.01 - 40 M ρ) ϕ^{n}_{i, j} + (10 M ρ - 0.99) ϕ^{n}_{i + 1, j} + (10 M ρ + 1) ϕ^{n}_{i - 1, j} + (10 M ρ +) ϕ_{i, j - 1}^{n} - 10 ϕ_{i, j}^{n + 1}

(3.18)

Take

i = 1, 2, 3, 4, 5

, ………

n = 0

j = 1

We get the algebraic equations

\begin{array}{l} - 20 M ρ ϕ^{o}_{2, 2} + (1 - 50 M ρ) ϕ^{0}_{1, 2} - 20 M ρ ϕ^{0}_{0, 2} = (10 M ρ - 0.99) ϕ^{0}_{2, 1} + (10.01 - 40 M ρ) ϕ^{0}_{1, 1} + (10 M ρ + 1) ϕ^{0}_{0, 1} + (10 M ρ + 1) ϕ^{0}_{1, 0} 0 - 10 ϕ^{1}_{1, 1} \\ - 20 M ρ ϕ^{o}_{32} + (1 - 50 M ρ) ϕ^{0}_{22} - 20 M ρ ϕ^{0}_{1, 2} = (10 M ρ - 0.99) ϕ^{0}_{3, 1} + (10.01 - 40 M ρ) ϕ^{0}_{2, 1} + (10 M ρ + 1) ϕ^{0}_{1, 1} + (10 M ρ + 1) ϕ^{0}_{2, 0} 0 - 10 ϕ^{1}_{2, 1} \\ - 20 M ρ ϕ^{o}_{4, 2} + (1 - 50 M ρ) ϕ^{0}_{3, 2} - 20 M ρ ϕ^{0}_{2, 2} = (10 M ρ - 0.99) ϕ^{0}_{4, 1} + (10.01 - 40 M ρ) ϕ^{0}_{3, 1} + (10 M ρ + 1) ϕ^{0}_{2, 1} + (10 M ρ + 1) ϕ^{0}_{3, 0} 0 - 10 ϕ^{1}_{3, 1} \\ - 20 M ρ ϕ^{o}_{5, 2} + (1 - 50 M ρ) ϕ^{0}_{4, 2} - 20 M ρ ϕ^{0}_{3, 2} = (10 M ρ - 0.99) ϕ^{0}_{5, 1} + (10.01 - 40 M ρ) ϕ^{0}_{4, 1} + (10 M ρ + 1) ϕ^{0}_{3, 1} + (10 M ρ + 1) ϕ^{0}_{4, 0} 0 - 10 ϕ^{1}_{4, 1} \\ - 20 M ρ ϕ^{o}_{6, 2} + (1 - 50 M ρ) ϕ^{0}_{5, 2} - 20 M ρ ϕ^{0}_{4, 2} = (10 M ρ - 0.99) ϕ^{0}_{6, 1} + (10.01 - 40 M ρ) ϕ^{0}_{5, 1} + (10 M ρ + 1) ϕ^{0}_{4, 1} + (10 M ρ + 1) ϕ^{0}_{5, 0} 0 - 10 ϕ^{1}_{5, 1} \end{array}

(3.19)

Taking the initial

ϕ (x, y, 0) = ϕ (x, 1, 0) = 1

, and boundary conditions,

ϕ (0, y, t) = ϕ (x, 0, t) = 1

, the above system of algebraic equation becomes

[\begin{matrix} (- 50 M ρ) & - 20 M ρ & 0 & 0 & 0 \\ - 20 M ρ & (- 50 M ρ) & - 20 M ρ & 0 & 0 \\ 0 & - 20 M ρ & (- 50 M ρ) & - 20 M ρ & 0 \\ 0 & 0 & - 20 M ρ & (- 50 M ρ) & - 20 M ρ \\ 0 & 0 & 0 & - 20 M ρ & (- 50 M ρ) \end{matrix}] [\begin{matrix} ϕ_{1, 2}^{0} \\ ϕ_{2, 2}^{0} \\ ϕ_{3, 2}^{0} \\ ϕ_{4, 2}^{0} \\ ϕ_{5, 2}^{0} \end{matrix}] = [\begin{matrix} 10 M ρ - 0.98 \\ - 10 M ρ + 1.02 \\ - 10 M ρ + 1.02 \\ - 10 M ρ + 1.02 \\ - 10 M ρ + 1.02 \end{matrix}]

(3.20)

Solving the above matrix equation, we get the solutions for varying values of

ρ

and M as shown in Table 1 and 2 in Section 4.

3.9. Reaction-diffusion Equation

We take into account a more involved reaction-diffusion model, the so-called Cahn-Hilliard model. This is the concept proposed by Chaplain and Anderson (2003), which emphasizes the mutually beneficial relationships between ECM and cancer tumour and metastatic capacities of cancer cells. We explore more elaborate extensions of the model and its various forms, including additional terms related to tumour cell proliferation, ECM renewal, and alternative functions representing tumour cell MDE generation. Our research aims to uncover how cancer cells are able to create and secrete the MDE, which in turn degrades the ECM and initiates migration of the cells towards healthy regions of the tissue. The model is given as;

\frac{\partial u}{\partial t} = D (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}}) + α u - β u

(3.21)

Additionally, to obtain the following simulations we use the following dimensionless parameter values:

D

α

and

β

3.10.1. Discretization of Reaction-diffusion Equation

For simplicity let D =1

\frac{M_{i, j}^{n + 1} - M_{i, j}^{n}}{Δ t} = \frac{M_{i + 1, j}^{n} - 2 M_{i, j}^{n} + M_{i - 1, j}^{n}}{{(Δ x)}^{2}} + \frac{M_{i, j + 1}^{n} - 2 M_{i, j}^{n} + M_{i, j - 1}^{n}}{{(Δ y)}^{2}} (α - β) M_{i, j}^{n}

(3.22)

Multiply by

Δ t

let

Δ x

Δ y

\frac{Δ t}{Δ x} = \frac{0.01}{0.1} = \frac{1}{0.1}, \frac{Δ t}{{(Δ x)}^{2}} = \frac{0.01}{{(0.1)}^{2}} = \frac{1}{100} = 0.01

M_{i, j + 1} - M_{i, j} = 0.01 M_{i + 1, j} - 0.02 M_{i, j} + 0.01 M_{i - 1, j} + 0.01 M_{i, j + 1} - 0.02 M_{i, j} + 0.01 M_{i, j + 1} + 0.01 (α - β) M_{i - 1, j}

(3.23)

0.01 M_{i + 1, j}^{n} + (1 - 0.04 + α - β) M_{i, j}^{n} + 0.01 M_{i - 1, j}^{n} = (1 - 0.01) M_{i, j + 1}^{n} - 0.01 M_{i, j - 1}^{n}

(3.24)

Simplifying, (3.24) gives

0.01 M_{i + 1, j}^{n} + (0.96 + α - β) M_{i, j}^{n} + 0.01 M_{i - 1, j}^{n} = 0.99 M_{i, j + 1}^{n} - 0.01 M_{i, j - 1}^{n}

(3.25)

In order for us to solve the system (3.25), we impose some initial conditions to get the algebraic equations below

\begin{array}{l} 0.01 M_{2, 1}^{0} + (0.96 + α - β) M_{1, 1}^{0} + 0.01 M_{0, 1}^{0} = 0.99 M M_{1, 2}^{0} - 0.01 M_{1, 0}^{0} \\ 0.01 M_{3, 1}^{0} + (0.96 + α - β) M_{2, 1}^{0} + 0.01 M_{1, 1}^{0} = 0.99 M M_{2, 2}^{0} - 0.01 M_{2, 0}^{0} \\ 0.01 M_{4, 1}^{0} + (0.96 + α - β) M_{3, 1}^{0} + 0.01 M_{2, 1}^{0} = 0.99 M M_{3, 2}^{0} - 0.01 M_{3, 0}^{0} \\ 0.01 M_{5, 1}^{0} + (0.96 + α - β) M_{4, 1}^{0} + 0.01 M_{3, 1}^{0} = 0.99 M M_{4, 2}^{0} - 0.01 M_{4, 0}^{0} \\ 0.01 M_{6, 1}^{0} + (0.96 + α - β) M_{5, 1}^{0} + 0.01 M_{4, 1}^{0} = 0.99 M M_{5, 2}^{0} - 0.01 M_{5, 0}^{0} \end{array}

(3.26)

Taking the initial

M (x, y, 0) = 10

and boundary conditions,

M (x, 0, 0) = 10

, the above system of algebraic equation becomes

[\begin{matrix} (0.96 + α - β) & 0.01 & 0 & 0 & 0 \\ 0.01 & (0.96 + α - β) & 0.01 & 0 & 0 \\ 0 & 0.01 & (0.96 + α - β) & 0.01 & 0 \\ 0 & 0 & 0.01 & (0.96 + α - β) & 0.01 \\ 0 & 0 & 0 & 0.01 & (0.96 + α - β) \end{matrix}] [\begin{matrix} M_{1, 1}^{0} \\ M_{2, 1}^{0} \\ M_{3, 1}^{0} \\ M_{4, 1}^{0} \\ M_{5, 1}^{0} \end{matrix}] = [\begin{matrix} 8.895 \\ 9.895 \\ 9.895 \\ 9.895 \\ 9.895 \end{matrix}]

(3.27)

Solving the above matrix equation, we get the solutions for varying values of

α

and

β

as shown in Table 4 and 5 in chapter 4 respectively.

3.11. Non-Dimensional Parameters

The study of dimensions is at the heart of dimensional analysis, a mathematical method used to provide answers to a wide range of engineering problems. In physics, there are two types of quantities: dimensional and non-dimensional, and any physical phenomenon can be described by an equation that describes the relationship between them. By integrating dimensional variables to generate non-dimensional parameters, dimensional analysis aids in establishing a systematic organization of the variables in the physical connection. Analyzing the complex fluid flow phenomenon is made possible by using dimensional analysis to plan model tests and display experimental results in a systematic manner. The following parameters play a significant role in our study.

3.9.1. Prandtl Number, Pr

The Prandtl number is a property of the fluid named after Ludwig Prandtl, a German mathematician associated with the concept of boundary layer theory. It is the ratio of viscous force to thermal force. It actually relates the relative thickness of the hydrodynamic and thermal boundary layers. Thus,

\Pr = \frac{v i s c o u s}{t h e r m a l} \frac{f o r c e}{f o r c e} = \frac{μ C p}{κ} . δ t

(3.28)

3.9.2. Reynolds Number, Re

Similar flow patterns can be predicted using this dimensionless number. It is the ratio of the inertial force to the viscous force. Mathematically written as:

Re = \frac{i n e r t i a}{v i s c o u s} \frac{f o r c e}{f o r c e} = \frac{ρ U_{0} H}{v}

(3.29)

The Reynolds number quantifies how much the viscous effect dominates the inertial effect. When the system's Re is low, the viscous force dominates, and viscosity has a significant impact on the flow field. If Re is big, however, the inertia force predominates, and the effect of viscosity is significant only in the localized region close to the boundary or anywhere else with a significant velocity gradient. The Reynolds number can be in a wide variety of values depending on the nature of the flow. Laminar flow is defined as a flow where the Reynolds number is less than 2000. Above 3500 indicates a turbulent flow, whereas values between 2000 and 3500 indicate a transitional flow. In the realm of fluid dynamics, it can be used to solve every problem imaginable. Reynolds number is applicable for closed surface flows as well as for free surface flows.

3.9.3. Peclet Number (Pe)

Jean-Claude Eugène Peclet of France is honoured with a dimensionless number that bears his name: the Peclet number. It is the ratio of the rate at which a physical quantity is advected by the flow to the rate at which the same quantity diffuses (either in matter or heat) under the influence of an adequate gradient. When discussing the transfer of species or mass, the Peclet number is calculated by multiplying the Reynolds number by the Schmidt number. The thermal Peclet number can be thought of as the product of the Reynolds and Prandtl numbers as they apply to thermal fluids. Fluid flows that involve both momentum and mass diffusion convection processes are defined by their Schmidt number, which is based on the Reynolds number to describe the flow regime.

P e_{L} = \frac{Advective transport rate}{mass diffusion rate} = {Re}_{L} S c = \frac{u}{D / L}

(3.30)

3.9.4. Schmidt number (Sc)

The Schmidt number, a dimensionless number named after German engineer Ernst Heinrich Wilhelm Schmidt (1892–1975). How quickly momentum may be transmitted from one location to another is measured by the kinematic viscosity, and this value is compared to the mass diffusivity to determine the Schmidt number (how easily mass can be transferred from one place to another). It's a way to describe flows of fluids where both momentum and mass diffusion convection mechanisms are active. The equations describing the mass momentum transfer are shown below.

S c = \frac{D m}{M d m}

(3.31)

D m

is the Molecular Diffusivity of Momentum and

M d m

is the Molecular Diffusivity of Mass. Molecular diffusivity of momentum is the spread of momentum between particles (atoms or molecules) of matter, often in the fluid state and Molecular diffusivity of mass is the spread of momentum between particles (atoms or molecules) of matter, often in the fluid state.

4. Results and Discussions

The simulation results given focus on the effects of the mobility, M, tumour density,

ρ

,and Peclet number, Pe, Extracellular Matrix,

α

and Matrix degradation Enzymes,

β

on tumour growth rate, tumour cells concentration and tumour invasion respectively.

4.1. Effects of Cell Density on Tumour Growth Rate

We solve equation (3.25) using MATLAB and obtain the results of the effects of

ρ

on tumour growth rate when we hold M constant and get results as shown in Table 4.1 below

The above results are presented in Figure 4.1 below

Table 4.1. Value of tumour growth rates for varying density.

Figure 4.1. Graph of tumour growth rate against length of tumour spread from source at varying tumour density.

Figure 4.1 shows how tumour density impacts tumour growth. Figure 4 shows that tumour growth slows with increasing tumour density and accelerates in the opposite direction as tumours spread from their original location. Tumours formed in high-density cells are suppressed through the application of compressive pressures, leading to smaller tumours than those grown in low-density cells. As a result, an increase in cell density typically results in a smaller tumour.

4.2. Effects of mobility on tumour growth rate

We solve equation (3.26) using MATLAB and get the results of the effects of M on tumour growth rate as shown in Table 4.2 below.

The above results are presented in Figure 4.2 below

Table 4.2. Value of tumour growth rates for Mobility.

Figure 4.2. Tumour growth rate against length of tumour spread from source at varying mobility.

Figure 4.2 shows how the movement of cancer cells affects the formation of tumours. Tumour growth rates are observed to rise alongside cell motility. As mobility decreases, the distance a tumour can travel from its initial site increases. Tumour and normal host cell motility play a role in tumour metastasis at multiple stages, including tumour cell motility during basement membrane breakdown, tumour cell motility during escape from the primary tumour, tumour cell motility during migration to blood and lymphatic vessels, and tumour cell motility during intravasation, extravasation, and metastasis to distant organs.

4.3. Effects of P’eclet Number on Tumour Concentration

We solve equation (3.14) using MATLAB and get the results of the effects of P’eclet number on tumour concentration as shown in Table 4.3 below

Table 4.3. Value of tumour concentration for various P’eclet number.

The above results are presented in Figure 4.3 below

Figure 4.3. Tumour concentration against length of tumour spread from source at varying P’eclet number.

The effect of P’eclet number on tumour concentration can be observed from Figure 6. Initially, an increase in P’eclet number leads to increases in the tumour growth rate. The effect of Pe on concentration is seen to have been very significant in two-phase flows in tumour cells. From the first phase, it is clearly seen that, with increasing values of Pe, the peak of tumour concentration increases from 0 to 1.5 and the profile of concentration becomes flatter. The peak values of tumour concentration are 1.5 to 2.0. With increasing values of Pe beyond 2.0, the initial tumour cells concentration decreases for a fixed radius from 2.0 to 4.0. From Figure 4.3 it is worth noting that the Pe governs the initial concentration of tumour cells in the human body system. For this reason, the concentration of the tumour cells reduces for higher values of Pe. The peak of the mean concentration decreases with the increase in the value of Pe and this effect is very significant in small cell radius. This study may be applied for investigation of the transportation process of drugs or plasma proteins in blood flow through small arteries.

4.4. Effects of ECM on Tumour Invasion

We solve equation (3.27) using MATLAB and get the results of the effects of ECM on tumour invasion rate as shown in Table 4.4 one below

Table 4.4. Value of tumour invasion for various ECM.

The above results are presented in Figure 4.4 below

Figure 4.4. Tumour invasion against length of tumour spread from source at varying ECM number.

Figure 4.4 displays the impact of extracellular matrix on tumour invasion. Rise in ECM leads to decreases in the tumour invasion rate. The ECM binds soluble substances, such as growth factors and other ECM-associated proteins. Skeletal muscle development from embryonic stage to senescence involves cell surface receptors interacting with ECM components and ECM-bound substances to regulate cell adhesion and cell signaling, consequently regulating activities as varied as proliferation, differentiation, migration, and apoptosis. Due to an increase in ECM in the cell milieu of aged muscle, satellite cells lose their capacity to differentiate into myogenic lineages. Collagen released by satellite cells can keep them dormant, and studies have shown that ECM protein components contribute to the myogenesis process of skeletal muscle progenitor cells. Hence, it is safe to say that ECM plays a crucial role in preventing tumour growth and sustaining normal cancer cell activity.

4.5. Effects of MDE on Tumour Invasion Rate

We solve equation (3.27) using MATLAB and get the results of the effects of MDE on tumour invasion rate as shown in Table 4.5 below

Table 4.5. Value of tumour invasion rate for various MDE.

The above results are presented in Figure 4.5 below

Figure 4.5. Tumour invasion rate against length of tumour spread from source at varying MDE.

Figure 4.5 shows how MDE affects the pace of tumour invasion. When DME levels rise, tumours invade at a faster rate. During DME, the invading tissue's extracellular matrix scaffolds are degraded in whole while a new, tumour-derived extracellular matrix is simultaneously formed to support the tumour mass's unchecked growth. As a result, it is safe to say that MDE plays a crucial role in promoting tumour growth by degrading cancer cells' physiological function and fostering the formation of skeletal muscle.

4.6. Equilibrium point between of ECM and MDEs in tumour invasion

We solve equation (3.27) using MATLAB and get the results of the point of equilibrium between ECM = 0.05 and MDE = 0.05 on tumour invasion rate as shown in Table 4.6 below.

Table 4.6. Value of Equilibrium point between of ECM and MDEs.

The above results are presented in Figure 4.6 below

Figure 4.6. Tumour invasion rate against length of tumour spread from source at equilibrium point between ECM and MDE.

Figure 4.6 demonstrates that the amount of extracellular matrix decreases with tumour size and distance from the tumour's origin. This is because the cells are stuck together, whereas the graph of matrix degrading enzymes shows a direct correlation between tumour length and tumour depth. Collagen, enzymes, and glycoproteins are just a few examples of the macromolecules and minerals that make up the extracellular matrix, which provides structural and biochemical support to the cells in their immediate vicinity. The integrins on the cells of the tissues it is connecting are the point of attachment for the extracellular matrix. Proteins known as integrins bind to the cytoskeleton of cells. Molecules of this kind are exudates from neighbouring cells. Yet, tumour cells' own matrix degrading enzymes (MDEs) break down the extracellular matrix, freeing cancer cells to spread across the body. As the point where the ECM and the MDEs meet, equilibrium is established.

5. Summary and Conclusions

A numerical study was performed to analyze the effects of the mobility M, tumour density,

μ

and Peclet number, Pe, Extracellular Matrix,

α

and Matrix degradation Enzymes,

β

on tumour growth rate tumour cells concentration and tumour invasion, respectively.

The following outcomes can be written from the study.

An increase in the mobility leads to an increase in the tumour growth rate.
Surface heat transfer tends to decrease with an increase in tumour growth rate
An increase in Peclet number leads to a decrease in tumour concentration.
An increase in ECM leads to an increase in tumour growth rate.
An increase in MDEs causes a decrease in tumour growth rate.
An increase in Peclet number leads to a decrease in tumour growth rate.

6. Recommendations

From this study, the following areas arise for further analysis and development:

An extension of this study to factor in a varying sources of tumour, ST and dead cells, SD.
An extension of this study to incorporate diffusion coefficient tumour growth rate

Funding Declaration:

The authors received no financial support for the research authorship and/or publication of this article.

Acknowledgement

The authors acknowledge Kisii University (Department of Mathematics and Actuarial Sciences) for affiliation and support in academic resources.

References

Akhtar, A., Majid, H., Abdul, G., Zafar, A., Kottakkaran, S., Alharthi M. and Wasim J. (2021). Numerical Simulations and Analysis for Mathematical Model of Avascular Tumour Growth Using Gompertz Growth Rate Function. Published by Elsevier BV. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/ licenses/by-nc-nd/4.0/).
Beerenwinkel, N., Roland F., Schwarz, Gerstung M., Markowetz F. (2015). Cancer Evolution: Mathematical models and computational inference, Systematic Biology, Volume 64, Issue 1, 1 January 2015, e25. [CrossRef]
Brassel, M. and Bretin, E. (2011). A modified phase field approximation for mean curvature flow with conservation of the volume, Math. Meth. Appl. Sci., 34, 1157–1180. [CrossRef]
Elaine, L. et al., (2009). Mathematical Oncology. Multiparameter Computational Modelling of Tumour Invasion. Available online: www.aacrjournals.org. [CrossRef]
Grego´rio, A., Fonseca, N., Moura, V., Lacierda, M., Figueredo, P., Simões, S., et al. (2016). Inoculated Cell Density as a Determinant Factor of the Growth Dynamics and Metastatic Efficiency of a Breast Cancer Murine Model. PLoS ONE 11(11): e0165817. [CrossRef]
Gross, J.F. and Popel, A.S. (1979). Mathematical models of transport phenomena in normal and neoplastic tissue, Tumour Blood Circulation. CRC Press, Peterson edition; (Google Scholar).
Hasitha, N., Weerasinghe, Pamela, M., Burrage, [...]and Dan V., Nicolau (Journal of Oncology). Published online 2019 oct 31.
Johnson, K.E., Howard. M. W, Strasser M.K., Lima EABF, Huang S, et al. (2019) Cancer cell population growth kinetics at low densities deviate from the exponential growth model and suggest an Allee effect. PLoS Biol 17(8). [CrossRef]
Katira, P., Roger, T., Bonnecaze and Muhammad H. Zaman., (2013). Modelling the mechanics of cancer: effect of changes in cellular and extra-cellular mechanical properties. Frontier in Oncology. [CrossRef]
Sharif, G.M. &Wellstein A. (2015). Cell density regulates cancer metastasis via the Hippo pathway. [CrossRef]
Socolofsky, S. A. and Jirka.G.H. (2002). Environmental Fluid Mechanics. Part 1, 2^nd Edition, Institute for Hydrodynamics, University of Karlsruhe, Germany.
Uthamacumaran, A. (2020). Cancer: A turbulence problem. Volume 22, Issue 12, Pages 759-769.Neoplasia. [CrossRef]
WHO (2020). Global Cancer Observatory, [2020, December 18]. Cancer data for cancer action. (http//gco.iarc.fr/today/home). GLOBOCAN_GCO [Twitter post] Retrieved from: http//twitter.com/today home.
Xuelian, H., Zhang L.Xuezhao, L… (2019). Toll-like receptor 2 and Toll-like receptor regulation of cancer cell stemness mediated by cell death-induced high-mobility group box 1. Biomedicine 40, 135–150. [CrossRef]
Zhang, X., Song, G., Liang, G. and Tian, T (2022). Mathematical Modelling and Analysis of Tumour Chemotherapy. Symmetry 2022, 14(4), 704. [CrossRef]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Mathematical Modeling of Dynamics of Cancer Invasion in Human Body Tissues

Abstract

1. Introduction

2. Background

2.1. Governing Equations and Theories

2.1.1. Continuity Equation

2.1.2. Momentum Equation

2.1.3. Concentration Equation

2.1.4. Energy Equation

2.1.5. Cahn–Hilliard Equation

3. Mathematical Formulation

3.4. Methods of Solution

3.4.1. Central Difference Scheme

3.4.2. Crank Nicolson Scheme

3.5. Dimensionalizing Tumour Concentration Equation

3.6. Discretization of Concentration Equation

3.7. Cahn-Hilliard Model

3.8. Discretization of Cahn-Hilliard Model

3.9. Reaction-diffusion Equation

3.10.1. Discretization of Reaction-diffusion Equation

3.11. Non-Dimensional Parameters

3.9.1. Prandtl Number, Pr

3.9.2. Reynolds Number, Re

3.9.3. Peclet Number (Pe)

3.9.4. Schmidt number (Sc)

4. Results and Discussions

4.1. Effects of Cell Density on Tumour Growth Rate

4.2. Effects of mobility on tumour growth rate

4.3. Effects of P’eclet Number on Tumour Concentration

4.4. Effects of ECM on Tumour Invasion

4.5. Effects of MDE on Tumour Invasion Rate

4.6. Equilibrium point between of ECM and MDEs in tumour invasion

5. Summary and Conclusions

6. Recommendations

Funding Declaration:

Acknowledgement

References

MDPI Initiatives

Important Links

Subscribe