# x t br or the

ρ(x, t) = √

or the curvature function

Our numerical experiments show that both the mesh density functions (19) and (20) produce reasonable adaptive meshes for solving the system (17) when the motility values Dc are between O(10−3) and O(10−2). However, for smaller motility

values, for example Dc = O(10−4), even though the cancer concentrations do not experience much variations at the early stage of the invasion, the solutions of the system (17) become quickly unstable due to the sharp structures that are developed as time evolves. This indicates that the magnitudes of the spatial gradients of the cancer concentrations become extremely large for this Diphenylterazine (DTZ) value of Dc . Therefore, for this size of motility value, the traditional arc-length mesh density function (19) becomes extremely large. In this case, the equidistribution of the arc-length mesh density function

(19) will result in mesh over concentration in these regions of large gradients and can potentially lead to the breakdown of the computation. Therefore, the mesh density function must be carefully chosen in order to avoid any mesh overlapping. Here, we propose the mesh density function to be defined as
ρ(x, t) =
(
1 + α

where α is a user defined parameter and

We find the mesh density function (21) to work well for all our computations using α = 5.

3.3. Solution procedure

Finally, we present the solution procedure employed for solving the cancer cell invasion model (18a)–(18e). We use an alternate solution procedure to solve the system (17) for the physical solution w and (8) for the steady solution and then taking its gradient to obtain the adaptive mesh. The steps of the alternate solution procedure are summarized as follows:

Given the adaptive mesh xn and the solution wn at time level n, then the adaptive mesh and the physical solution are solved alternately as follows:

4. Numerical details

We use the method of lines (MOL) for the numerical solution of the system (17). We first approximate all the spatial

derivatives in (17) in terms of the computational domain variables. The resulting system of ordinary differential equations (ODEs) is integrated using a suitable ODE solver. In what follows we give the details of the spatial discretization for the case of one spatial dimension, the discretization of the 2D case follows similarly.

D
ψξ

xξ

xξ

xξ

where D represents the diffusion coefficients Dc , Du, Dp and Dm and ψ represents the solution components c, u, p and m.
The chemotaxis terms in (18a) are approximated using central upwind discretizations described as follows:
xξ
[cK ]ξ j
= 2
xj+1
− xj−1
.
(24)

where the gradients (cξ )j are computed using the nonlinear flux limiter scheme

The minmod function is defined as

if

minj
φj
,

otherwise.

The convection terms in (17) are discretized using the upwind scheme (see [35] for more details)
x wξ

x

wj

xj

xj

Notice that the boundary points are fixed, i.e. x˙ = 0 for j = 1 and j = N. We would like to point out that the upwind scheme (30) helps to stabilize the solution of the physical problem as mentioned in [35].

dw
w

where f (wj, t) denotes the spatial discretization of (18a)–(18e). We would like to remark that a small time step ∆t is required when using the alternate solution procedure described in Section 3. This is required so that the adaptive mesh does not change significantly within the time interval [t, t + ∆t]. In our computations, we find that choosing ∆t = 0.01 the system (31) can be integrated efficiently using MATLAB ode113 solver, a variable time step Adams–Bashforth–Moulton method [36].