<figure><img src="image_1.png" /></figure>

HAL
open science

# A BDF B-spline Method for a Nonlocal Tumor Growth Model

Abderrahman Bouhamidi, El Harraki Imad, Yassine Melouani

<figure><img src="image_2.png" /></figure>

## To cite this version:

Abderrahman Bouhamidi, El Harraki Imad, Yassine Melouani. A BDF B-spline Method for a Nonlocal Tumor Growth Model. Computers & Mathematics with Applications, In press. hal-04994394

HAL Id: hal-04994394

https://hal.science/hal-04994394v1

Submitted on 17 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A BDF B-spline Method for a Nonlocal Tumor Growth Model

Bouhamidi Abderrahman$^a$, El Harraki Imad$^b$, Melouani Yassine$^{a,b}$

$^a$Université du littoral côte d'opale, 50 rue Ferdinand Buisson, CS 80699, Calais, 62100, France

$^b$Ecole nationale supérieure des mines de Rabat, Av. Hadj Ahmed Cherkaoui - B.P. : 753, Rabat, 10000, Morocco

**Abstract**

This paper presents a model for tumor growth using nonlocal velocity. We establish some results on the existence and uniqueness of the solution for a nonlocal tumor growth model. Many experiences show that tumor spheroid can be invariant by rotation and can guard the shape of the spheroid during the growth process in some particular cases. Here, we assume that the multiple components of the system are invariant by rotation. Then, we use the Backward Differentiation Formula (BDF) spline to solve the nonlocal system. To illustrate the efficiency of the proposed method, we performed numerical tests that simulate a tumor growth scenario. Such techniques may be used to provide informations on practical applications of the model.

**Keywords:** B-spline Collocation Method, Backward Differentiation Formula, Nonlocal Tumor Growth Model, Partial Differential Equations

**Introduction**

Mathematical modeling of tumor growth is a promising area of research that has a potential to improve our understanding of cancer and develop more effective treatments. By simulating the behavior of cancer cells and their interactions with their environment, researchers can gain insights into the factors that contribute to tumor growth and optimize treatment strategies. In fact, one of the primary motivations for studying tumor growth through mathematical modeling is to gain a deeper understanding of the underlying biological processes that drive cancer, such as genetic mutations, immune system response, and blood vessel formation. This understanding can then be used to identify potential targets for new cancer treatments. In recent years, an important increase in tumor growth models has been studied, see for instance [1, 2, 3, 4, 5, 6] and references therein. However, a better understanding of this phenomenon is still needed due to the diversity of complex factors that play a role in cancer invasion and metastasis.

Recent studies have shown that cancer cells can communicate through a variety of mechanisms, including direct cell-to-cell contact, secretion of signaling molecules, and extracellular vesicles [7, 8]. This communication allows cells to share information about their environment and coordinate their behavior to promote tumor growth and survival. Furthermore, there is increasing evidence that collective cell migration is common during the invasion and metastasis of malignant tumors. This behavior has been observed and investigated in many biological areas, for example, predator spray systems, fishing, flocks of birds, and insect colonies [9].

Most classical approaches to modeling tumor growth using balance equations or diffusive partial differential equations do not take these collective patterns into account [1, 10]. However, experience over the past several years has shown that nonlocal balance laws present a very good way to model the behavior of collective and nonlocal movement of the population. In [11], a nonlocal predator-prey model has been studied and in [12, 13], a class of nonlocal pedestrian models has been experienced. Nonlocal modeling has been used also in economy; we refer here to some nonlocal supply-chain models [10, 14].

In this paper, we are going to extend the model studied in [4], to a nonlocal balance laws model. Some similar models have been investigated in [15, 16, 17, 18, 19, 20] and references therein, where in [15, 16, 17, 18] a nonlocal diffusion operator has been used, while in [19, 20] the nonlocality has been laid on the advection term using a nonlocal velocity.

This paper is presented as following. In Section 1, we introduce a simplified nonlocal model for tumor growth, laying the groundwork for the subsequent analysis. Section 2 is dedicated to establishing the existence and uniqueness of
the solution, leveraging general results from [21] and [22]. In Section 3, we focus on a radial model derived from the initial multidimensional model, examining its specific characteristics and implications. Section 4 devoted to the construction of a numerical scheme using the Cubic B-splines method for spatial discretization and the Backward Differentiation Formula (BDF) for temporal integration. Finally, in Section 5, we perform accuracy tests to illustrate the efficiency of the proposed method and perform a numerical tests that simulate tumor growth scenarios in the absence of medical treatment. This allows to provide informations on practical applications of the model.

# 1. A nonlocal continuum model

Mathematical modeling of solid tumor growth using continuum models has been extensively studied, see [23, 3, 24, 25], and references therein. This type of modeling is compatible with solid and secondary tumors, also for multicellular spheroids.

The general form of the equation describing the tumor behavior is as follows:

$$ \frac{\partial P}{\partial t}(t, x) + \operatorname{div}(v(t, x)P(t, x)) = \text{Production} - \text{Death}, \qquad (1) $$

where $(t, x) \in [0, T] \times \mathbb{R}^d, T > 0, d \ge 2, v$ is the velocity vector, and $P$ is the density of proliferation cells. In the practical tumor growth modeling the dimension $d$ is usually 1, 2 or 3. The velocity vector is commonly determined by using Darcy's law under the consideration of tumor cells as incompressible flow. The downside of this equation is that it does not take into account the collective behavior of tumor cells, cell-to-cell adhesion, and other types of regrouping [26]. Thus, we use a nonlocal type of modeling, taking into account the adhesion in one hand, between tumor cells, and in the other hand between healthy cells. We consider the following nonlocal equation of tumor growth:

$$ \frac{\partial P}{\partial t}(t, x) + \operatorname{div}[V[\vec{\alpha}_p, \gamma_p, P](t, x)P(t, x)] = H(M(t, x))P(t, x) - a_1\lambda(t, x)M(t, x)P(t, x), \qquad (2) $$

where $\vec{\alpha}_p$ is the directional vector of the velocity, $\gamma_p$ is the kernel function of the proliferation, $M$ represents the concentration of oxygen and nutriment in the tissue, $\lambda$ is the density of the drug, while $a_1$ is the resistance of the drug by tumor cells and the function $H$ represents the growth factor. Let us remark that if $a_1 = 0$, the tumor is resistant to this type of medical treatment. Here we assume that the proliferation cells accept the treatment, and so we have $a_1 > 0$.

Here, the flux is assumed to be nonlocal, so the velocity vector $V$ is expressed as follows:

$$ V[\vec{\alpha}_p, \gamma_p, P](t, x) = \vec{\alpha}_p(\gamma_p * P)(t, x) = \vec{\alpha}_p(t, x) \iint_{\mathbb{R}^d} \gamma_p(x - y)P(t, y)dy. \qquad (3) $$

As usual, see [19], the kernel function $\gamma_p$ is chosen as a Gaussien type function. The production rate is controlled by the growth factor $H$, and by the medical treatment $\lambda$.

As in [27, 4], we assume that the growth function $H$ is given as follows:

$$ H(M) = \kappa \frac{1 + \tanh(\delta(M - M_{threshold}))}{2} $$

This function indicates that if the supplement concentration is greater than the hypoxia threshold $M_{threshold}$, the tumor grows by a factor of $\kappa$, while if this amount is below the hypoxia value, growth is suspended. The positive constant $\delta$ controls the growth of the hyperbolic tangent function.

We add here a second equation describing the movement of healthy cells under medical treatment, as follows:

$$ \frac{\partial S}{\partial t}(t, x) + \operatorname{div}[V[\vec{\alpha}_p, \gamma_p, P](t, x)P(t, x)] = -a_2\lambda(t, x)M(t, x)S(t, x), \qquad (4) $$

where $S$ represents the density of healthy cells, $a_2$ is a constant that denotes the secondary effects of medical treatment, $\vec{\alpha}_s$ is a radial directional vector, $\gamma_s$ is a Gaussian kernel. The nonlocal velocity for healthy cells is presented as follows:

$$ V[\vec{\alpha}_s, \gamma_s, S](t, x) = \vec{\alpha}_s(t, x) (\gamma_s * S). \qquad (5) $$

The concentration of nutrients and oxygen is a crucial detail in tumor growth. The increase in this concentration is controlled by several factors. One of the main factors is angiogenesis, which is a fundamental biological process
through which new blood vessels are formed from preexisting ones. In the context of tumor development, angiogenesis is often stimulated by signaling molecules and growth factors produced by cancer cells themselves. This process leads to the creation of an extensive network of blood vessels that infiltrate the tumor mass. This process can be omitted here under the condition that tumor cells have not yet reached the vascular stage. However, healthy cells are capable of producing and stimulating new resources of oxygen and nutrients when needed. However, excessive consumption of oxygen and nutrients by cancer cells causes a lack of supplements in the area.

Finally, we add a third equation that describes the concentration of oxygen and nutrients as following:

$$ \frac{\partial M}{\partial t}(t, x) - D\Delta M(t, x) = M_s S(t, x)(1 - M(t, x)) - \eta M(t, x)P(t, x), \quad (6) $$

where $M_s$ is the rate at which healthy cells can create new blood vessels, $\eta$ is the rate of consumption of oxygen and nutrients by cancer cells, while $D$ denotes the diffusion of the vascularisation in the tissue.

We assume that the domain of study is large enough such that there is no interaction between tumor cells and the boundary. In light of that and without loss of generality, we assume that the domain of our study is the open disk of center zero and radius $R > 0$, denoted as

$$ \Omega_R = \{x \in \mathbb{R}^d | \|x\|_2 < R\}, $$

where $\| \cdot \|_2$ is the usual euclidean norm in $\mathbb{R}^d$, and let $\partial\Omega_R$ be the boundary of $\Omega_R$. As there is no interaction of tumor cells with the boundary, we have

$$ P = 0 \quad \text{in} \quad \partial\Omega_R. $$

In the other hand we assume that the density of healthy cells $S$ and the concentration of oxygen and nutrients $M$ remain constant throughout the boundary; consequently we have

$$ S = S_R, \quad M = M_R \quad \text{in} \quad \partial\Omega_R. $$

With $S_R$ and $M_R$ are two fixed positive values.

We restrict the convolution terms $\gamma_p * P$ and $\gamma_s * S$ inside the computational domain, therefor we replace 3 and 5 by

$$ V[\vec{\alpha}_p, \gamma_p, P](t, x) = \vec{\alpha}_p(t, x) \int_{\mathbb{R}^d} 1_{\Omega_R} \gamma_p(x-y) P(t, y) dy = \vec{\alpha}_p(t, x) \int_{\Omega_R} \gamma_p(x-y) P(t, y) dy \quad \text{and} \quad (7) $$

$$ V[\vec{\alpha}_s, \gamma_s, S](t, x) = \vec{\alpha}_s(t, x) \int_{\mathbb{R}^d} 1_{\Omega_R} \gamma_s(x-y) S(t, y) dy = \vec{\alpha}_s(t, x) \int_{\Omega_R} \gamma_s(x-y) S(t, y) dy \quad (8) $$

To ensure that the growth process remains in $\Omega_R$ and there is no interference with the boundary, let $R' > 0$ with $R' < R$, there exists a function $m: \mathbb{R}^d \to \mathbb{R}$ that is infinitely differentiable, satisfying:

$$ m(x) = \begin{cases} 1 & \text{if } \|x\|_2 \le R', \\ 0 & \text{if } \|x\|_2 \ge R. \end{cases} $$

We multiply the right-hand side of the equations of proliferative cells and healthy ones by the function $m$. We regroup equations (2), (4), and (6) with their boundary conditions as follows:

$$ \left\{
\begin{array}{ll}
\displaystyle \frac{\partial P}{\partial t}(t, x) + \text{div} \left[ V[\vec{\alpha}_p, \gamma_p, P](t, x) P(t, x) \right] = m(x) \left( H(M(t, x)) P(t, x) - a_1 \lambda(t) M(t, x) P(t, x) \right), & (t, x) \in [0, T] \times \Omega_R \\[2ex]
\displaystyle \frac{\partial S}{\partial t}(t, x) + \text{div} \left[ V[\vec{\alpha}_p, \gamma_p, P](t, x) P(t, x) \right] = -a_2 m(x) \lambda(t) M(t, x) S(t, x), & (t, x) \in [0, T] \times \Omega_R \\[2ex]
\displaystyle \frac{\partial M}{\partial t}(t, x) - D \Delta M(t, x) = M_s S(t, x)(1 - M(t, x)) - \eta M(t, x) P(t, x), & (t, x) \in [0, T] \times \Omega_R \\[2ex]
\displaystyle P(t, x) = 0, \quad S(t, x) = S_R, \quad M(t, x) = M_R, & (t, x) \in [0, T] \times \partial \Omega_R \\[2ex]
\displaystyle P(0, x) = P_0(x), \quad S(0, x) = S_0(x), \quad M(0, \cdot) = M_0(x), & x \in \Omega_R.
\end{array}
\right.
\quad (9) $$
Let $U = (P, S, M)$, we rewrite the system of equations (9) into the following form:

$$ egin{cases} 
\frac{\partial U}{\partial t}(t, x) + \mathcal{A}(U)(t, x) = G(U)(t, x), & (t, x) \in [0, T] \times \Omega_R \\ 
U(t, x) = U_R, & (t, x) \in [0, T] \times \partial \Omega_R \\ 
U(0, x) = U_0(x). & x \in \Omega_R 
\end{cases} 
\quad (10) $$

with

$$ \mathcal{A}(U)(t, x) = \left[ \operatorname{div} \left( \vec{\alpha}_p(t, x) \left( \int_{\Omega_R} \gamma_p(x-y) P(t, y) dy \right) P(t, x) \right), \operatorname{div} \left( \vec{\alpha}_s(t, x) \left( \int_{\Omega_R} \gamma_s(x-y) S(t, y) dy \right) S(t, x) \right), -D\Delta M \right], $$

$$ G(U) = \begin{pmatrix} g_1(U) \\ g_2(U) \\ g_3(U) \end{pmatrix} = \begin{pmatrix} mP(H(M) - a_1\lambda M) \\ -ma_2\lambda MS \\ m_sS(1-M) - \eta MP \end{pmatrix}, \quad U_R = (0, S_R, M_R) \quad \text{and} \quad U_0 = (P_0, S_0, M_0). $$

## 2. Existence and uniqueness of solution

Our aim in this section is to prove the existence and uniqueness of solutions for the system outlined in equation (9). We begin by examining the local dynamics through the isolation of the convolution term and leveraging established principles from the theory of semilinear evolution equations. Following this, we integrate insights from the theory of nonlocal balance equations, as elaborated in [21] and [22], to affirm the existence and uniqueness of the solutions for the system mentioned in equation (2).

Before proceeding with the existence and uniqueness proofs, we need to define the following spaces.

Let $L^\infty(\Omega_R)$ be the space of essentially bounded measurable functions on $\Omega_R$, equipped with the norm:

$$ \|f\|_{L^\infty(\Omega_R)} = \ess\sup_{x \in \Omega_R} |f(x)|. $$

Let $C(\Omega_R)$ be the space of continuous functions on $\Omega_R$, with the uniform norm:

$$ \|f\|_{C(\Omega_R)} = \sup_{x \in \Omega_R} |f(x)|. $$

Let $C_b^1(\Omega_R)$ be the space of continuously differentiable functions on $\Omega_R$ with bounded derivatives, normed by:

$$ \|f\|_{C_b^1(\Omega_R)} = \|f\|_{C(\Omega_R)} + \|\nabla f\|_{C(\Omega_R)}. $$

For a Banach space $(X, \|\cdot\|_X)$, let $C([0, T]; X)$ denote the space of continuous functions from $[0, T]$ to $X$, with the norm:

$$ \|f\|_{C([0,T];X)} = \sup_{t \in [0,T]} \|f(t)\|_X. $$

Let $L^1([0, T]; X)$ be the space of Bochner integrable functions from $[0, T]$ to $X$, with the norm:

$$ \|f\|_{L^1([0,T];X)} = \int_0^T \|f(t)\|_X dt. $$

Let $W^{2,1}(\Omega_R)$ be the Sobolev space defined as:

$$ W^{2,1}(\Omega_R) = \{u \in L^1(\Omega_R) : D^\alpha u \in L^1(\Omega_R) \text{ for all } |\alpha| \le 2\} $$

### 2.1. Existence and uniqueness of the solution of the local system

Let $w_p$ and $w_s$ be two fixed functions belonging to $C([0, T], C_b^1(\mathbb{R}^d))$, and let $w := (w_p, w_s)$. We consider the associated local system to (10) and written as follows:

$$ egin{cases} 
\frac{\partial U_w}{\partial t}(t, x) + \mathcal{A}_w(U_w)(t, x) = G(U_w)(t, x), & (t, x) \in [0, T] \times \Omega_R \\ 
U_w(t, x) = U_R, & (t, x) \in [0, T] \times \partial \Omega_R \\ 
U(0, x) = U_0(x). & x \in \Omega_R 
\end{cases} 
\quad (11) $$
with $A_w(U) = [\text{div}(\vec{\alpha}_p(t, x)w_p(t, x)P_w(t, x)), \text{div}(\vec{\alpha}_s(t, x)w_s(t, x)S_w(t, x)), -D\Delta M]$, and

$$G(U_w) = \begin{pmatrix} g_1(U_w) \\ g_2(U_w) \\ g_3(U_w) \end{pmatrix} = \begin{pmatrix} mP_w(H(M_w) - a_1\lambda M_w) \\ -ma_2\lambda M_w S_w \\ M_s S_w(1 - M_w) - \eta M_w P_w \end{pmatrix}$$

and we have $U_R = (0, S_R, M_R)$ and $U_0 = (P_0, S_0, M_0)$.

To ensure that the solution remains positive and bounded within the domain, we make the following assumptions regarding initial conditions and velocities.

**Assumptions.** We assume that:

$$(A1) \quad U_0 \in C^1(\Omega_R)^3 := C^1(\Omega_R) \times C^1(\Omega_R) \times C^1(\Omega_R),$$

$$(A2) \quad U_0 = (P_0, S_0, M_0) \ge 0, \text{ which means } P_0 \ge 0, S_0 \ge 0, \text{ and } M_0 \ge 0,$$

$$(A3) \quad U_0 - U_R \text{ is compactly supported in } \Omega_R,$$

$$(A4) \quad \vec{\alpha}_k \in C^2([0, T] \times \Omega_R; \mathbb{R}^d) \quad k = p, s.$$

$$(A5) \quad \vec{\alpha}_s \text{ is compactly supported in } [0, T] \times \Omega_R,$$

$$(A6) \quad \lambda \in C([0, T]).$$

The Banach space $L^\infty(\Omega_R)^3 := L^\infty(\Omega_R) \times L^\infty(\Omega_R) \times L^\infty(\Omega_R)$ is equipped here with the norm:

$$\|f\|_{L^\infty(\Omega_R)^3} = \|f_1(t)\|_{L^\infty(\Omega_R)} + \|f_2(t)\|_{L^\infty(\Omega_R)} + \|f_3(t)\|_{L^\infty(\Omega_R)} \quad \text{for all } f = (f_1, f_2, f_3) \in L^\infty(\Omega_R)^3.$$

**Theorem 1.** Let the assumptions (A1)-(A6) hold. Then there exists a maximal time $T_{max}$ such that equation (11) has a unique positive solution in $C([0, T], L^\infty(\Omega_R)^3)$. Furthermore, $U(t, \cdot)_w - U_R$ is a compactly supported function in $\Omega_R$ for all $t \le T_{max}$.

*Proof.* To establish the existence and uniqueness of (11), we apply the characteristic method. First, we ensure that there is no interaction between the tumor and the boundary by selecting $T$ so that the characteristics are appropriately defined within $\Omega_R$. This is guaranteed by the fact that $U_0 - U_R$ and the velocity vector of healthy cells have a compact support in $\Omega_R$. Let $0 < R_0 < R$ be such that $\text{supp}(U_0) \cup \text{supp}(\vec{\alpha}_p) \subset \Omega_{R_0}$,

and we let

$$R_i = R_0 + \frac{i(R - R_0)}{3} \quad \text{for } i = 1, 2, 3. \qquad (12)$$

In order to prevent the characteristic curves from exiting $\Omega_R$, we proceed under the assumption that $T \le T_{max}$ is sufficiently small such that

$$T < T_1^* := \frac{R - R_0}{3 \max\{\|\vec{\alpha}_p w_p\|_{C([0, T_{max}], L^\infty(\Omega_R))}, \|\vec{\alpha}_s w_s\|_{C([0, T_{max}], L^\infty(\Omega_R))}\}} \quad (13)$$

with $T_{max}$ is a large fixed value. The characteristics are well defined for all $t, s \in [0, T]$ and $x \in \Omega_{R_2}$ as follows

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial X[t, x]_{w_p}(s)}{\partial s} = \vec{\alpha}_p(s, X[t, x]_{w_p}(s))w_p(s, X[t, x]_{w_p}(s)), \\ X[t, x]_{w_p}(t) = x. \end{array} \right. \qquad (14)$$

And

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial X[t, x]_{w_s}(s)}{\partial s} = \vec{\alpha}_s(s, X[t, x]_{w_s}(s))w_p(s, X[t, x]_{w_s}(s)), \\ X[t, x]_{w_s}(t) = x. \end{array} \right. \qquad (15)$$

In fact, let $x \in \Omega_{R_2}$. By using the condition on the final time $T$ we have for $k = p, s$:

$$\begin{align*} \|X[t, x]_{w_k}(s)\|_2 &\le \|x\|_2 + \|X[t, x]_{w_k}(s) - x\|_2 \\ &\le R_2 + \int_t^s \|\vec{\alpha}_k(s, X[t, x]_{w_k}(s))w_p(s, X[t, x]_{w_k}(s))\|_2 ds \\ &\le R_0 + \frac{2}{3}(R - R_0) + T \|\vec{\alpha}_k w_k\|_{C([0, T], L^\infty(\Omega_R))} \\ &< R_0 + \frac{2}{3}(R - R_0) + \frac{1}{3}(R - R_0) = R. \end{align*}$$
Using the semi-explicit form for the local balance law in [22] we have for $(t, x) \in [0, T] \times \Omega_{R_2}$

$$P_w(t, x) = P_0(X[t, x]_{w_p}(0)) \det(DX[t, x]_{w_p}(0)) + \int_0^T \det(DX[t, x]_{w_p}(s)) g_1(U_w) (s, X[t, x]_{w_p}(s)) ds$$

$$S_w(t, x) = S_0(X[t, x]_{w_s}(0)) \det(DX[t, x]_{w_s}(0)) + \int_0^T \det(DX[t, x]_{w_s}(s')) g_2(U_w) (s', X[t, x]_{w_s}(s')) ds'$$

with

$$\det(DX[t, x]_{w_k}(s)) = \exp\left(\int_t^s \operatorname{div}(\alpha_k(s, X[t, x]_{w_k}(s))w_p(s, X[t, x]_{w_k}(s))) ds\right) \quad \text{for } k = p, s. \quad (16)$$

Now, letting $t, s \in [0, T]$, and letting $x \in \Omega_{R_2} \setminus \Omega_{R_1}$, we have then

$$\|X[t, x]_{w_k}(s)\|_2 \ge \|x\|_2 - \|X[t, x]_{w_k}(s) - x\|_2, \quad k = p, s \\ > \left|R_1 - \frac{1}{3}(R - R_0)\right| = R_0.$$

Using the assumption that initial functions are compactly supported in $\Omega_{R_0}$, we deduce that $P_0(X[t, x]_{w_p}(0)) = 0$ and $S_0(X[t, x]_{w_s}(0)) = S_R$.

Furthermore, the function $m$ ensures that the right-hand side is compactly supported in $\Omega_R$, by choosing $R' \le R_0$ and using the assumption of the compactness of the support of $\alpha_s$, we get $P_w(t, x) = 0$ and $S_w(t, x) = S_R$ for all $(t, x) \in [0, T] \times \Omega_{R_2} \setminus \Omega_{R_1}$. Using equation (11), we can smoothly extend $P$ and $S$ to be constant in $\Omega_R \setminus \Omega_{R_2}$.

For the equation of nutrients $M$, it is a semi-linear parabolic equation, therefore, by using the maximum principle we deduce that $U_w - U_R$ is compactly supported in $\Omega_R$. Now we prove the local existence and uniqueness for Equation (11). To do this, we employ a classical approach for semilinear evolutionary equations using Banach's fixed-point theorem. This approach is demonstrated in works such as [28, 29].

In the first step, we show that the right hand side is locally lipschitzien. In fact, let $L > 0$, for $U_w = (P_w, S_w, M_w)$ and $V_w = (P'_w, S'_w, M'_w)$ in $C([0, T], L^\infty(\Omega_R)^3)$ such that $\|U_w\|_{C([0,T],L^\infty(\Omega_R)^3)}, \|V_w\|_{C([0,T],L^\infty(\Omega_R)^3)} \le L$, we have
$$\|G(U_w(t))-G(V_w(t))(t)\|_{L^\infty(\Omega_R)^3} = \|g_1(U_w(t))-g_1(V_w(t))\|_{\infty}+\|g_2(U_w(t))-g_2(V_w(t))\|_{\infty}+\|g_3(U_w(t))-g_3(V_w(t))\|_{\infty}.$$

By using the fact that $\|H(M)\|_{C([0,T],L^\infty(\Omega_R))} \le \kappa$ together with following inequalities:

$$\|g_1(U_w(t)) - g_1(V_w(t))\|_{\infty} \le \kappa \|P_w(t) - P'_w(t)\|_{\infty} + a_1 L (\|P_w(t) - P'_w(t)\|_{\infty} + \|M_w(t) - M'_w(t)\|_{\infty}), \\ \|g_2(U_w(t)) - g_2(V_w(t))\|_{\infty} \le a_2 L (\|S_w(t) - S'_w(t)\|_{\infty} + \|M_w(t) - M'_w(t)\|_{\infty}) + M_s M_R \|S_w(t) - S'_w(t)\|_{\infty},$$

and

$$\|g_2(U_w(t))-g_2(V_w(t))\|_{\infty} \le M_s L (\|S_w(t) - S'_w(t)\|_{\infty} + \|M_w(t) - M'_w(t)\|_{\infty}) + \eta L (\|P_w(t) - P'_w(t)\|_{\infty} + \|M_w(t) - M'_w(t)\|_{\infty}), \\ \text{we get}$$

$$\|G(U_w(t)) - G(V_w(t))\|_{L^\infty(\Omega_R)^3} \le C(L)\|U_w(t) - V_w(t)\|_{\infty}, \quad (17)$$

where $C(L)$ is a positive function depending on the constant $L$. Finally, by taking the maximum over time $t \in [0, T]$, we get the result.

For the next step, we consider the following constants:

$$C_k = \|\alpha_k\|_{\infty} \|\nabla w_k\|_{\infty} + \|\operatorname{div}(\alpha_k)\|_{\infty} \|w_k\|_{\infty} \quad \text{and} \quad C_v = e^{T \max\{C_p, C_s\}} \ge 1. \quad (18)$$

Let $K = 2LC_v$. We denote by

$$E = \{f \in C([0, T], L^\infty(\Omega_R)^3) \mid \|f\|_{C([0,T],L^\infty(\Omega_R)^3)} \le K\}, \quad (19)$$

the closed ball of radius $K$ (closed convex) subset of the Banach space $C([0, T], L^\infty(\Omega_R)^3)$.

Let $U_w \in E$, we consider the following functional

$$\Phi(U_w)(t) = \begin{pmatrix} P_0(X[t, x]_{w_p}(0)) \det(DX[t, x]_{w_p}(0)) + \int_0^T \det(DX[t, x]_{w_p}(s)) g_1(U_w) (s, X[t, x]_{w_p}(s)) ds \\ S_0(X[t, x]_{w_s}(0)) \det(DX[t, x]_{w_s}(0)) + \int_0^T \det(DX[t, x]_{w_s}(s')) g_2(U_w) (s', X[t, x]_{w_s}(s')) ds' \\ T(t)M_0 + \int_0^T T(t-s)g_2(U_w)(s)ds \end{pmatrix},$$
where $t \to \mathcal{T}(t)$ is the semigroup of the heat operator $A = -D\Delta$, with the domain:

$$D(A) = \{u \in W^{2,1}(\Omega_R) : Au \in L^\infty(\Omega_R), u|_{\partial\Omega_R} = M_R\}.$$

The semigroup $t \to \mathcal{T}(t)$ satisfies that:

$$\|\mathcal{T}(t)M_0\|_{C([0,T],L^\infty(\Omega_R))} \le \|M_0\|_{L^\infty(\Omega_R)} \quad \text{for all } t \in [0,T], \qquad (20)$$

see [30] for more details. First, let us show that $\Phi(E) \subseteq E$. Let $U_w \in E$, we have

$$\begin{aligned} \|\Phi(U_w)(t)\|_{\infty} = & \left\| P_0(X[t, x]_{w_p}(0)) \det(DX[t, x]_{w_p}(0)) + \int_0^T \det(DX[t, x]_{w_p}(s)) g_1(U_w)(s, X[t, x]_{w_p}(s)) ds \right\|_{\infty} \\ & + \left\| S_0(X[t, x]_{w_s}(0)) \det(DX[t, x]_{w_s}(0)) + \int_0^T \det(DX[t, x]_{w_s}(s')) g_2(U_w)(s', X[t, x]_{w_s}(s')) ds' \right\|_{\infty} \\ & + \left\| \mathcal{T}(t)M_0 + \int_0^T \mathcal{T}(t-s)g_3(U_w)(s)ds \right\|_{\infty}. \end{aligned}$$

By using (16),(18),(20) and estimates from [22], we obtain

$$\|\det(DX[t, x]_{w_k}(s))\|_{C([0,T], L^\infty(\Omega_R))} \le C_v \quad \text{for } k = p, s. \qquad (21)$$

Then, it follows that

$$\|\Phi(U_w)(t)\|_{\infty} \le C_v (\|P_0\|_{\infty} + \|S_0\|_{\infty}) + \|M_0\|_{\infty} + T[C_v (\|g_1(U_w)(t)\|_{\infty} + \|g_2(U_w)(t)\|_{\infty}) + \|g_3(U_w)(t)\|_{\infty}]. \quad (22)$$

Using the lipschitz property of $G$ and $G(0) = 0$, we have

$$\begin{aligned} \|\Phi(U_w)(t)\|_{\infty} &\le C_v [(\|P_0\|_{\infty} + \|S_0\|_{\infty} + \|M_0\|_{\infty}) + TC(L) (\|g_1(U_w)(t)\|_{\infty} + \|g_2(U_w)(t)\|_{\infty} + \|g_3(U_w)(t)\|_{\infty})] \\ &\le C_v L(1 + TC(L)). \end{aligned}$$

Then, by choosing

$$T \le T_2^* := \min \left\{ T_1^*, \frac{1}{2C(L)} \right\}, \qquad (23)$$

and taking the maximum over time, we deduce that $\Phi(E) \subseteq E$.

Now, let $U_w, U'_w \in E$, we have

$$\|\Phi(U_w) - \Phi(U'_w)\|_E \le TC_v \|U_w - U'_w\|_E \le \frac{1}{2} \|(U_w - U'_w)\|_E.$$

Therefore, $\Phi$ is a contraction in $E$ with Lipschitz constant $1/2$, and so $\Phi$ has a fixed point $U_w \in E$, which ensures the existence and uniqueness of a solution of (11). $\square$

## 2.2. Existence and uniqueness of the nonlocal system

Now we turn back into the nonlocal system (2), we use Banach's fixed theorem once again as in [22], therefore, we announce the main theorem

**Theorem 2.** Let the assumptions (A1)-(A6) hold, and assume that $\gamma_p, \gamma_s \in C_b^1(\Omega_R)$. Then the system of equations (4) has a unique solution $U = (P, S, M)$ in $C([0, T], L^\infty(\Omega_R)^3)$.

*Proof.* In the proof, we follow the lines of ([22]) adapted to the case of a system of coupled nonlocal equations. Let define the following constants:

$$N_p := \|\gamma_p\|_{C([0,T] \times \Omega_R)} \left( \|P_0\|_{L^1(\Omega_R)} + \sup_{w \in C([0,T], C_b^1(\Omega_R))^2} \left( \|g_1(U_w)\|_{L^1([0,T], L^1(\Omega_R))} \right) \right) \quad \text{and}$$

$$dN_p := \left\| \frac{d\gamma_p}{dx} \right\|_{C([0,T] \times \Omega_R)} \left( \|P_0\|_{L^1(\Omega_R)} + \sup_{w \in C([0,T], C_b^1(\Omega_R))^2} \left( \|g_1(U_w)\|_{L^1([0,T], L^1(\Omega_R))} \right) \right),$$
and

$$
\begin{align*}
N_s &:= \|\gamma_s\|_{C([0,T] \times \Omega_R)} \left( \|S_0\|_{L^1(\Omega_R)} + \sup_{w \in C([0,T], C_b^1(\Omega_R))^2} \left( \|g_2(U_w)\|_{L^1([0,T], L^1(\Omega_R))} \right)^2 \right), \\
dN_s &:= \left\| \frac{d\gamma_s}{dx} \right\|_{C([0,T] \times \Omega_R)} \left( \|S_0\|_{L^1(\Omega_R)} + \sup_{w \in C([0,T], C_b^1(\Omega_R))^2} \left( \|g_2(U_w)\|_{L^1([0,T], L^1(\Omega_R))} \right)^2 \right).
\end{align*}
$$

$$
\text{Where } C([0, T], C_b^1(\Omega_R))^2 := C([0, T], C_b^1(\Omega_R)) \times C([0, T], C_b^1(\Omega_R)).
$$

Remark. We have shown in Theorem 1 that by the compactness of initial data, the velocities, and the semilinear terms, we have the boundedness of the local solution, therefor $G(U_w) = (g_1(U_w), g_2(U_w), g_3(U_w))$ is bounded for every $w \in C([0, T], C_b^1(\Omega_R))^2$.

We set $N = \max\{N_p, N_s\}$ and $N_z = \max\{dN_p, dN_s\}$. Then we consider

$$
B = \left\{ f \in C([0, T], C_b^1(\Omega_R))^2 \mid \|f\|_{C([0,T], C(\Omega_R))^2} \le N, \quad \|\nabla f\|_{C([0,T], C(\Omega_R))^2} \le N_z \right\} \quad (24)
$$

It is known that $B$ is closed subset of a Banach space. Now we check that $\mathcal{F}$ is a contraction on $B$.
Let $w := (w_p, w_s) \in B$, we consider the following mapping $\mathcal{F}$ as follow

$$
\mathcal{F}(w)(t, x) = \begin{pmatrix} \mathcal{F}_1(w)(t, x) \\ \mathcal{F}_2(w)(t, x) \end{pmatrix} = \begin{pmatrix} \int_{X_{w_p}[t, \Omega_R](0)} \gamma_p(x - X_{w_p}[0, y](t)) P_0(y) dy + \int_{\Omega_R} \int_0^T \gamma_p(x - X_{w_p}[s, y](t)) g_1(U_w)(s, y) dyds \\ \int_{X_{w_s}[t, \Omega_R](0)} \gamma_s(x - X_{w_s}[0, y](t)) S_0(y) dy + \int_{\Omega_R} \int_0^T \gamma_s(x - X_{w_s}[s', y](t)) g_2(U_w)(s', y) dyds' \end{pmatrix}. \tag{25}
$$

Let us introduce the following notations:

$$
\Lambda_1(P,g) := \|P\|_{L^1(\Omega_R)} + \sup_{w \in C([0,T], C_b^1(\Omega_R))^2} \left( \|g(U_w)\|_{L^1([0,T], L^1(\Omega_R))} \right),
$$

$$
\Lambda_{\infty}(P,g) := \|P\|_{L^{\infty}(\Omega_R)} + \sup_{w \in C([0,T], C_b^1(\Omega_R))^2} \left( \|g(U_w)\|_{L^1([0,T], L^{\infty}(\Omega_R))} \right),
$$

$$
C_{\infty}(P, g, \gamma) := 1 + 8RC_v \|\gamma\|_{C([0,t] \times \Omega_R)} \Lambda_{\infty}(P, g).
$$

Following the same demonstration as done in [22], we can show that $\mathcal{F}(B) \subset B$ and we have for $w, w' \in B$ the
following estimates:

$$
|\mathcal{F}_1(w)(t,x) - \mathcal{F}_1(w')(t,x)| \le \left\| X_{w_p}[t, \cdot](*) - X_{w'_p}[t, \cdot](*) \right\|_{C([0,t] \times \Omega_R)} \|\gamma_p(t, \cdot)\|_{C_b^1(\Omega_R)} \Lambda_1(P_0, g_1) C_\infty(P_0, g_1, \gamma_p).
$$

As

$$
\left\| X_{w_p}[t, \cdot](*) - X_{w'_p}[t, \cdot](*) \right\|_{C([0,t] \times \Omega_R)} \le T \|w_p - w'_p\|_{C([0,t]; C(\Omega_R))} \exp \left( T \left( N_z + \|D\vec{\alpha}_p\|_{L^\infty([0,t] \times \Omega_R)} \right) \right),
$$

we get

$$
\begin{equation}
\begin{split}
|\mathcal{F}_1(w)(t,x) - \mathcal{F}_1(w')(t,x)| \le{}& T \|w_p - w'_p\|_{C([0,t];C(\Omega_R))} \exp\left(T(N_z + \|D\vec{\alpha}_p\|_{L^\infty([0,t]\Omega_R)})\right) \|\gamma_p(t, \cdot)\|_{C_b^1(\Omega_R)} \Lambda_1(P_0, g_1) C_\infty(P_0, g_1, \gamma_p).
\end{split}
\tag{26}
\end{equation}
$$

Using the same estimates, we get similarly

$$
\begin{equation}
\begin{split}
|\mathcal{F}_2(w)(t,x) - \mathcal{F}_2(w')(t,x)| \le{}& T \|w_s - w'_s\|_{C([0,t];C(\Omega_R))} \exp\left(T(N_z + \|D\vec{\alpha}_s\|_{L^\infty([0,t]\Omega_R)})\right) \|\gamma_s(t, \cdot)\|_{C_b^1(\Omega_R)} \Lambda_1(S_0, g_2) C_\infty(S_0, g_2, \gamma_s).
\end{split}
\tag{27}
\end{equation}
$$

By summing the two estimates (26) and (27), we deduce that $\mathcal{F}$ is Lipschitzien, then, by choosing $T$ small enough, we have

$$
\|\mathcal{F}(w) - \mathcal{F}(w')\|_{C([0,T], C(\Omega_R))^2} < \frac{1}{2} \|w - w'\|_{C([0,T], C(\Omega_R))^2} \quad (28)
$$
Therefore, $\mathcal{F}$ is a contraction on $B$, then using Banach's fixed point theorem, there exists a unique fixed point $w \in B$ such that $\mathcal{F}(w) = w$ for $t \le T^*$. Furthermore, [31, Theorem 3.24] ensures that the existence and the uniqueness of the solution hold also for any final time $T > 0$, hence, in our case we have the existence and uniqueness for all $t \le T_2^*$ given in (23). $\square$

### 3. Radial model reformulation

The experience shows that the tumor spheroid can be invariant by rotation and can guard the shape of the spheroid during the growth process in some cases. So, we assume that the tumor remains in spheroid form for all $t \le T$, and we use the assumption of rotational invariance to write the model in radial coordinates, which allows us to obtain a simplified expression of the model.

Let $f: \mathbb{R}^d \to \mathbb{R}$ a radial function, by definition, it exists a function which we note $\tilde{f}$ satisfies

$$f(x) = \tilde{f}(\|x\|) = \tilde{f}(r) \quad \text{for all } x \in \mathbb{R}^d, \text{ with } r = \|x\|.$$

Under the assumption of invariance by rotation, $P$ is a radial function, we assume also that $S, M, \alpha_k, \gamma_k$, and $m$ are radial for $k = p, s$.

The next proposition is crucial for the characterisation of a radial expression for the nonlocal term.

**Proposition 1** (The convolution of two radial functions). Let $f$ and $g$ be two radial functions, defined from $\mathbb{R}^d$ to $\mathbb{R}$ such that $f * g$ is well defined. Then the convolution product $f * g$ is also a radial function. Furthermore, we have for all $x \in \mathbb{R}^d$:

$$f * g(x) = \tilde{f} * \tilde{g}(r) = \tilde{f} * \tilde{g}(r), \qquad (29)$$

where $\tilde{f} * \tilde{g}$ is defined as

$$\tilde{f} * \tilde{g}(r) = \frac{2\pi^{\frac{d-1}{2}}}{\Gamma(\frac{d-1}{2})} \int_{0}^{+\infty} \left[ \int_{0}^{\pi} \tilde{f} \left( \sqrt{s^2 + r^2 + 2rs \cos(\theta)} \right) s^{d-1} \sin^{d-2}(\theta) d\theta \right] \tilde{g}(s) ds, \qquad (30)$$

where $\Gamma$ denotes the classical Gamma function.

*Proof.* Let $x_1$ and $x_2$ two $\mathbb{R}^d$ vectors, such that $\|x_1\| = \|x_2\|$. Then it exists an orthogonal endomorphism $A$ such that $Ax_2 = x_1$, then we have

$$f * g(x_1) = \int_{\mathbb{R}^d} f(x_1 - y)g(y)dy = \int_{\mathbb{R}^d} f(Ax_2 - y)g(y)dy.$$

With the substitution $y = Az$, using the fact that $|\det(A)| = 1$, we get

$$f * g(x_1) = \int_{\mathbb{R}^d} f(Ax_2 - Az)g(Az)dz.$$

So we have

$$f * g(x_1) = \int_{\mathbb{R}^d} f(A(x_2 - z))g(Az)dz = \int_{\mathbb{R}^d} f(x_2 - z)g(z)dz = f * g(x_2).$$

Which shows that $f * g$ is a radial function. By using the fact that the convolution product $f * g$ is radial, we get

$$\begin{aligned} f * g(x) = f * g(\|x\|e_1) &= \int_{\mathbb{R}^d} f(\|x\|e_1 - y)g(y)dy = \int_{\mathbb{R}^d} \tilde{f}(\|\|x\|e_1 - y\|)\tilde{g}(\|y\|)dy \\ &= \int_{\mathbb{R}^d} \tilde{f}\left(\sqrt{(\|x\| - y_1)^2 + \|\tilde{y}\|_2^2}\right)\tilde{g}\left(\sqrt{y_1^2 + \|\tilde{y}\|_2^2}\right)dy, \end{aligned}$$

with $e_1 = (1, 0, \dots, 0) \in \mathbb{R}^d$ and $\tilde{y} = (y_2, y_3, \dots, y_d) \in \mathbb{R}^{d-1}$. Using the Fubini theorem, we get the following:

$$f * g(x) = \int_{\mathbb{R}} \left( \int_{\mathbb{R}^{d-1}} \tilde{f}\left(\sqrt{(\|x\| - y_1)^2 + \|\tilde{y}\|_2^2}\right) \tilde{g}\left(\sqrt{y_1^2 + \|\tilde{y}\|_2^2}\right) d\tilde{y} \right) dy_1.$$
We recall that for a radial function h on $\mathbb{R}^d$, we have the formula:

$$ \int_{\mathbb{R}^d} h(x)dx = \omega_d \int_0^{+\infty} \tilde{h}(r)r^{d-1}dr, $$

where $\omega_d = \frac{2\pi^{d/2}}{\Gamma(\frac{d}{2})}$ denotes the measure of the unit sphere in $\mathbb{R}^d$. Then, we may write

$$ f * g(x) = \omega_{d-1} \int_{\mathbb{R}} \left( \int_{0}^{+\infty} \tilde{f} \left( \sqrt{(\|x\| - y_1)^2 + \tilde{r}^2} \right) \tilde{g} \left( \sqrt{y_1^2 + \tilde{r}^2} \right) \tilde{r}^{d-2} d\tilde{r} \right) dy_1, $$

where $\tilde{r} = \|\tilde{y}\|$ is 2-norm in $\mathbb{R}^{d-1}$. Using the polar coordinates:

$$ \Psi: ]0, +\infty[ \times [ -\frac{\pi}{2}, \frac{\pi}{2} [ \rightarrow ]0, +\infty[ \times \mathbb{R} \\ (s, \varphi) \mapsto (\tilde{r}(s, \varphi) = s \cos(\varphi), y_1(s, \varphi) = s \sin(\varphi)), $$

we obtain

$$ \begin{aligned} f * g(x) &= \omega_{d-1} \int_{0}^{+\infty} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \tilde{f} \left( \sqrt{(\|x\| - s \sin(\varphi))^2 + s^2 \cos^2(\varphi)} \right) \tilde{g}(s) s^{d-2} \cos^{d-2}(\varphi) \, sd\varphi ds \\ &= \omega_{d-1} \int_{0}^{+\infty} \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}} \tilde{f} \left( \sqrt{\|x\|^2 - 2s\|x\|^2 \sin(\varphi) + s^2} \right) \tilde{g}(s) s^{d-1} \cos^{d-2}(\varphi) \, d\varphi ds \end{aligned} $$

And by using the following change of variable $\theta = \varphi + \frac{\pi}{2}$, we get the desired result. $\square$

Now we will transform the model (9) into radial coordinates, due to the fact that the solution is compactly supported, see Theorem 1, then we have

$$ \begin{aligned} \mathrm{div} (V[\vec{\alpha}_p, \gamma_p, P](t, x)P(t, x)) &= \mathrm{div} (\vec{\alpha}_p (\gamma_p * (1_{\Omega_R} P))(t, x)P(t, x)) = \mathrm{div} (\vec{\alpha}_p (\gamma_p * P)(t, x)P(t, x)) \\ &= \mathrm{div} (\vec{\alpha}_p(t, x) (\gamma_p * P)(t, x)) P(t, x) + (\vec{\alpha}_p(t, x) (\gamma_p * P)(t, x)) \cdot \nabla P(t, x) \\ &= \left[ \mathrm{div} (\vec{\alpha}_p(t, x)) (\gamma_p * P)(t, x) + \vec{\alpha}_p(t, x) \cdot \nabla \left( (\gamma_p * P)(t, x) \right) \right] P(t, x) + \\ &\quad + (\vec{\alpha}_p(t, x) (\gamma_p * P)(t, x)) \cdot \nabla P(t, x), \end{aligned} $$

The directional vector of the velocity $\vec{\alpha}_p$, is assumed to be a radial vector, so it can be written as follows

$$ \vec{\alpha}_p(t, x) = \widetilde{\vec{\alpha}}_p(t, r) \vec{e}_r, $$

with $\vec{e}_r = \frac{x}{\|x\|}$. So we have

$$ \mathrm{div} (\vec{\alpha}_p(t, x)) = \frac{d-1}{r} \widetilde{\alpha}_p(t, r) + \frac{\partial \widetilde{\alpha}_p(t, r)}{\partial r} = \frac{1}{r^{d-1}} \frac{\partial}{\partial r} (r^{d-1} \widetilde{\alpha}_p(t, r)), $$

and as

$$ \nabla ((\gamma_p * P)(t, x)) = \frac{\partial(\widetilde{\gamma_p * P})}{\partial r}(t, r) \vec{e}_r, \quad \text{and} \quad \nabla P(t, x) = \frac{\partial \tilde{P}}{\partial r}(t, r) \vec{e}_r, $$

we get the following:

$$ \vec{\alpha}_p(t, x) \cdot \nabla ((\gamma_p * P)(t, x)) = \widetilde{\alpha}_p(t, r) \frac{\partial(\widetilde{\gamma_p * P})}{\partial r}(t, r) = \widetilde{\alpha}_p(t, r) \frac{\partial \widetilde{\gamma_p}}{\partial r} * \widetilde{P}(t, r), $$

and

$$ \vec{\alpha}_p(t, x) (\gamma_p * P)(t, x) \cdot \nabla P(t, x) = \widetilde{\alpha}_p(t, r) (\widetilde{\gamma_p * P})(t, r) \frac{\partial \tilde{P}}{\partial r}(t, r). $$
Thus, we have

$$
\begin{align*}
\operatorname{div}(V[\vec{\alpha}_p, \gamma_p, P](t, x)P(t, x)) = & \left[ \frac{1}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \widetilde{\alpha}_p(t, r) \right) (\widetilde{\gamma}_p * \widetilde{P})(t, r) + \widetilde{\alpha}_p(t, r) \left( \frac{\partial \widetilde{\gamma}_p}{\partial r} * \widetilde{P} \right)(t, r) \right] \widetilde{P}(t, r) \\
& + \widetilde{\alpha}_p(t, r) (\widetilde{\gamma}_p * \widetilde{P})(t, r) \frac{\partial \widetilde{P}}{\partial r}(t, r) \\
& = \widetilde{\alpha}_p(t, r) \frac{\partial ((\widetilde{\gamma}_p * \widetilde{P}) \widetilde{P})}{\partial r}(t, r) + \frac{1}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \widetilde{\alpha}_p(t, r) \right) (\widetilde{\gamma}_p * \widetilde{P})(t, r) \widetilde{P}(t, r).
\end{align*}
$$

Finally we can write the equation for proliferation cells in the system (9) in the following form

$$
\frac{\partial \tilde{P}}{\partial t}(t, r) + \tilde{\alpha}_p(t, r) \frac{\partial ((\tilde{\gamma}_p * \tilde{P}) \tilde{P})}{\partial r}(t, r) = \left( \tilde{m}(r) H(\tilde{M}(t, r)) - a_1 \tilde{m}(r) \lambda(t) \tilde{M}(t, r) - \frac{1}{r^{d-1}} \frac{\partial}{\partial r} (r^{d-1} \tilde{\alpha}_p(t, r)) (\tilde{\gamma}_p * \tilde{P})(t, r) \right) \tilde{P}(t, r). \quad (31)
$$

In the same way, for the equation of healthy cells, we have

$$
\frac{\partial \tilde{S}}{\partial t}(t, r) + \tilde{\alpha}_s(t, r) \frac{\partial ((\tilde{\gamma}_s * \tilde{S}) \tilde{S})}{\partial r}(t, r) = \left( -a_2 \tilde{m}(r) \lambda(t) \tilde{M}(t, r) - \frac{1}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \tilde{\alpha}_s(t, r) \right) (\tilde{\gamma}_s * \tilde{S}) (t, r) \right) \tilde{S}(t, r). \quad (32)
$$

We recall that

$$
\Delta M(t, x) = \Delta \tilde{M}(t, \|x\|_2) = \frac{d-1}{r} \frac{\partial \tilde{M}}{\partial r}(t, r) + \frac{\partial^2 \tilde{M}}{\partial r^2}(t, r) = \frac{1}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \frac{\partial \tilde{M}}{\partial r}(t, r) \right). \quad (33)
$$

Finally, we write the system (9) in radial coordinates as follows

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial \tilde{P}}{\partial t}(t,r) + \tilde{\alpha}_p(t,r) \frac{\partial ((\tilde{\gamma}_p * \tilde{P}) \tilde{P})}{\partial r}(t,r) = \left( \tilde{m}(r) H(\tilde{M}(t,r)) - a_1 \tilde{m}(r) \lambda(t) \tilde{M}(t,r) - \frac{1}{r^{d-1}} \frac{\partial}{\partial r} (r^{d-1} \tilde{\alpha}_p(t,r)) (\tilde{\gamma}_p * \tilde{P})(t,r) \right) \tilde{P}(t,r), \\
\\
\displaystyle \frac{\partial \tilde{S}}{\partial t}(t,r) + \tilde{\alpha}_s(t,r) \frac{\partial ((\tilde{\gamma}_s * \tilde{S}) \tilde{S})}{\partial r}(t,r) = \left( -a_2 \tilde{m}(r) \lambda(t) \tilde{M}(t,r) - \frac{1}{r^{d-1}} \frac{\partial}{\partial r} (r^{d-1} \tilde{\alpha}_s(t,r)) (\tilde{\gamma}_s * \tilde{S})(t,r) \right) \tilde{S}(t,r), \quad (t,r) \in [0, T] \times [0, R] \\
\\
\displaystyle \frac{\partial \tilde{M}}{\partial t}(t,r) - \frac{D}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \frac{\partial \tilde{M}}{\partial r}(t,r) \right) = M_s \tilde{S}(t,r) (1 - \tilde{M}(t,r)) - \eta \tilde{M}(t,r) \tilde{P}(t,r), \quad (t,r) \in [0, T] \times [0, R], \\
\\
P(t,R) = 0, \quad S(t,R) = S_R, \quad M(t,R) = M_R, \qquad t \in [0, T], \\
\\
P(0,r) = P_0(r), \quad S(0,r) = S_0(r), \quad M(0,r) = M_0(r), \qquad r \in ]0, R].
\end{array}
\right.
\tag{34}
$$

We denote by $\tilde{U} = (\tilde{P}, \tilde{S}, \tilde{M})$ the exact solution of the radial problem (34). Then the system (34) may be written
in the following form:

$$
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial \tilde{U}}{\partial t}(t,r) + \tilde{\mathcal{A}}(\tilde{U})(t,r) = \tilde{G}(\tilde{U})(t,r), & (t,r) \in [0, T] \times [0, R], \\
\\
\qquad \tilde{U}(t,R) = \tilde{U}_R, & t \in [0, T], \\
\\
\qquad \tilde{U}(0,r) = \tilde{U}_0(r). & r \in ]0, R],
\end{array}
\right.
\tag{35}
$$

with

$$
\tilde{\mathcal{A}}(\tilde{U})(t,r) = \left[ \tilde{\alpha}_p(t,r) \frac{\partial ((\tilde{\gamma}_p * \tilde{P}) \tilde{P})}{\partial r}(t,r), \tilde{\alpha}_s(t,r) \frac{\partial ((\tilde{\gamma}_s * \tilde{S}) \tilde{S})}{\partial r}(t,r), -\frac{D}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \frac{\partial \tilde{M}}{\partial r}(t,r) \right) \right],
$$
$$
\tilde{G}(\tilde{U})(t,r) = (\tilde{g}_1(\tilde{U})(t,r), \tilde{g}_2(\tilde{U})(t,r), \tilde{g}_3(\tilde{U})(t,r)),
$$

where

$$
\begin{align*}
\tilde{g}_1(\tilde{U})(t,r) &= \begin{cases} \tilde{m}(r)H(\tilde{M}(t,r)) - a_1\tilde{m}(r)\lambda(t)\tilde{M}(t,r) - \frac{1}{r^{d-1}}\frac{\partial}{\partial r}(r^{d-1}\tilde{\alpha}_p(t,r))(\tilde{\gamma}_p * \tilde{P})(t,r) \tilde{P}(t,r), \\ -a_2\tilde{m}(r)\lambda(t)\tilde{M}(t,r) - \frac{1}{r^{d-1}}\frac{\partial}{\partial r}(r^{d-1}\tilde{\alpha}_s(t,r))(\tilde{\gamma}_s * \tilde{S})(t,r) \tilde{S}(t,r), \end{cases} \\
\tilde{g}_2(\tilde{U})(t,r) &= \begin{cases} \tilde{m}(r)H(\tilde{M}(t,r)) - a_1\tilde{m}(r)\lambda(t)\tilde{M}(t,r) - \frac{1}{r^{d-1}}\frac{\partial}{\partial r}(r^{d-1}\tilde{\alpha}_p(t,r))(\tilde{\gamma}_p * \tilde{P})(t,r) \tilde{P}(t,r), \\ -a_2\tilde{m}(r)\lambda(t)\tilde{M}(t,r) - \frac{1}{r^{d-1}}\frac{\partial}{\partial r}(r^{d-1}\tilde{\alpha}_s(t,r))(\tilde{\gamma}_s * \tilde{S})(t,r) \tilde{S}(t,r), \end{cases} \\
\tilde{g}_3(\tilde{U})(t,r) &= \frac{\partial \tilde{M}}{\partial t}(t,r) - \frac{D}{r^{d-1}}\frac{\partial}{\partial r}\left(r^{d-1}\frac{\partial \tilde{M}}{\partial r}(t,r)\right) = M_s \tilde{S}(t,r)(1-\tilde{M}(t,r)) - \eta \tilde{M}(t,r) \tilde{P}(t,r), \\
U_R &= (0, S_R, M_R) \text{ and } U_0 = (P_0, S_0, M_0).
\end{align*}
$$
## 4. The BDF B-spline method

In this section, we provide a description of the collocation method used to solve the radial system (35). The collocation method utilizes the cubic B-splines, which are widely favored in numerical analysis due to their inherent smoothness and flexibility in function approximation. B-splines are piecewise polynomial functions that ensure continuity up to the second derivative, making them particularly effective for interpolating data points or solving differential equations.

To begin with, we uniformly divide the spatial interval $[0, R]$ into $N + 7$ points. This division is crucial for applying the B-spline method as it specifies the nodes where the B-splines will be computed. The subdivision points are labeled as follows:

$$r_{-3} < r_{-2} < r_{-1} < 0 < r_{0} < r_{1} < \dots < r_{N} = R < r_{N+1} < r_{N+2} < r_{N+3},$$

where $[r_0, R]$ is our domain of study, with $r_0$ is a very small value, and the nodes $r_i$ are defined by the relation $r_i = ih$ for $i = -3, \dots, N + 3$ with the step size $h$ given by $h = \frac{R - r_0}{N}$. The points $r_{-3}, r_{-2}, r_{-1}, r_{N+1}, r_{N+2}, r_{N+3}$ are added to manage boundary conditions and ensure the spline approximation remains accurate and stable near the edges of the interval. In this context, the fundamental B-spline function refers to the cubic spline centered at the nodes $-2, -1, 0, 1, 2$, and its support is limited to the interval $[-2, 2]$. It can be written as follows:

$$B(r) = \begin{cases} 0 & \text{if } r < -2 \text{ or } r \ge 2, \\ \frac{1}{6}(2+r)^3 & \text{if } -2 \le r < -1, \\ \frac{1}{6}(4-6r^2-3r^3) & \text{if } -1 \le r < 0, \\ \frac{1}{6}(4-6r^2+3r^3) & \text{if } 0 \le r < 1, \\ \frac{1}{6}(2-r)^3 & \text{if } 1 \le r < 2, \end{cases} \quad (36)$$

and $B_i$ represent the well-known B-spline functions associated to the nodes $r_i$ defined for $i = -1, \dots, N + 1$ by

$$B_i(r) = B\left(\frac{r - r_i}{h}\right).$$

The function $B_i$ is supported by the interval $[r_{i-2}, r_{i+2}]$. Table 1 provides a summary of the values of $B_i$ and its derivatives at the points $r_i$.

<table><caption>**Table 1:** Values of B-splines and their derivatives at the points $r_i$.</caption><thead><tr><th>x</th><th>r<sub>i-2</sub></th><th>r<sub>i-1</sub></th><th>r<sub>i</sub></th><th>r<sub>i+1</sub></th><th>r<sub>i+2</sub></th></tr></thead><tbody><tr><td>B<sub>i</sub>(r)</td><td>0</td><td>1/6</td><td>4/6</td><td>1/6</td><td>0</td></tr><tr><td>B'<sub>i</sub>(r)</td><td>0</td><td>-1/2h</td><td>0</td><td>1/2h</td><td>0</td></tr><tr><td>B''<sub>i</sub>(r)</td><td>0</td><td>1/h<sup>2</sup></td><td>-2/h<sup>2</sup></td><td>1/h<sup>2</sup></td><td>0</td></tr></tbody></table>

The exact solution of (35) is approximated by a cubic B-spline given by the following expression:

$$\tilde{U}_h(t, r) = (\tilde{P}_h(t, r), \tilde{S}_h(t, r), \tilde{M}_h(t, r)) = \left( \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t) B_i(r), \sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) B_i(r), \sum_{i=-1}^{N+1} \alpha_i^{(3)}(t) B_i(r) \right)$$

where $\alpha_i^{(k)}(t)$ are time-dependent coefficients with unknown values for $k = 1, 2, 3$. Consider the following vector-
valued functions of dimensions $(N-1) \times 1$ and $3(N-1) \times 1$ respectively:

$$ \mathbb{B}(r) = \begin{pmatrix} B_1(r) \\ \vdots \\ B_{N-1}(r) \end{pmatrix} \quad \text{and} \quad \alpha(t) = \begin{bmatrix} \alpha^{(1)}(t) \\ \alpha^{(2)}(t) \\ \alpha^{(3)}(t) \end{bmatrix} = \begin{pmatrix} \alpha_1^{(1)}(t) \\ \vdots \\ \alpha_{N-1}^{(1)}(t) \\ \alpha_1^{(2)}(t) \\ \vdots \\ \alpha_{N-1}^{(2)}(t) \\ \alpha_1^{(3)}(t) \\ \vdots \\ \alpha_{N-1}^{(3)}(t) \end{pmatrix} \quad (37) $$

The function $U_h(t, x)$ has the following form

$$ \tilde{U}_h(t,r) = \begin{pmatrix} \tilde{P}_h(t,r) \\ \tilde{S}_h(t,r) \\ \tilde{M}_h(t,r) \end{pmatrix} = \begin{pmatrix} \alpha_{-1}^{(1)}(t)B_{-1}(r) + \alpha_{0}^{(1)}(t)B_{0}(r) + \mathbb{B}(r)^T \alpha^{(1)}(t) + \alpha_{N}^{(1)}(t)B_{N}(r) + \alpha_{N+1}^{(1)}(t)B_{N+1}(r) \\ \alpha_{-1}^{(2)}(t)B_{-1}(r) + \alpha_{0}^{(2)}(t)B_{0}(r) + \mathbb{B}(r)^T \alpha^{(2)}(t) + \alpha_{N}^{(2)}(t)B_{N}(r) + \alpha_{N+1}^{(2)}(t)B_{N+1}(r) \\ \alpha_{-1}^{(3)}(t)B_{-1}(r) + \alpha_{0}^{(3)}(t)B_{0}(r) + \mathbb{B}(r)^T \alpha^{(3)}(t) + \alpha_{N}^{(3)}(t)B_{N}(r) + \alpha_{N+1}^{(3)}(t)B_{N+1}(r) \end{pmatrix}, \quad (38) $$

where the notation $\mathbb{B}(r)^T$ stands for the transpose of the vector $\mathbb{B}(r)$. The collocation method involves substituting $\tilde{U}_h$ and its derivatives in (35) with the expression for $\tilde{U}_h$ given by (38). This substitution is followed by evaluating the resulting equation at the points $r_i$ for $i = 0, \dots, N$. Consequently, we obtain an ordinary differential equation that we can solve later. Beginning with the first equation for tumor cells $\tilde{P}_h$, we have:

$$ \frac{\partial \tilde{P}_h}{\partial t}(t,r) + \tilde{\alpha}_p(t,r) \frac{\partial \left( (\tilde{\gamma}_p * \tilde{P}_h) \tilde{P}_h \right)}{\partial r}(t,r) = \tilde{g}_1(\tilde{U}_h)(t,r), $$

which may be written in the form

$$ \frac{\partial \tilde{P}_h}{\partial t}(t,r) + (\tilde{\gamma}_p * \tilde{P}_h)(t,r) \left[ \frac{1}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \tilde{\alpha}_p(t,r) \right) \tilde{P}_h(t,r) + \tilde{\alpha}_p(t,r) \frac{\partial \tilde{P}_h}{\partial r}(t,r) \right] + \left( \frac{\partial \tilde{\gamma}_p}{\partial r} * \tilde{P}_h \right)(t,r) \tilde{P}_h(t,r) = (H(\tilde{M}_h(t,r)) - a_1 \lambda(t) \tilde{M}_h(t,r)) \tilde{P}_h(t,r), \quad (39) $$

We first start by expressing the nonlocal term, then

$$ (\tilde{\gamma}_p * \tilde{P}_h)(t,r) = \left( \tilde{\gamma}_p * \left( \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t) B_i \right) \right) (t,r) = \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t) (\tilde{\gamma}_p * B_i)(t) \quad (40) $$

We note for $i = -1, \dots, N+1$, and $k = p, n$:

$$ Z_i^k(r) := (\tilde{\gamma}_p * B_i)(r), \quad \text{and} \quad \bar{Z}_i^k(r) := \frac{\partial (\tilde{\gamma}_p * B_i)(r)}{\partial r} = \left( \frac{\partial \tilde{\gamma}_p}{\partial r} * B_i \right)(r). \quad (41) $$

Then by substituting in (40), we get

$$ \sum_{i=-1}^{N+1} \frac{d\alpha_i^{(1)}}{dt}(t)B_i(r) + \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t)Z_i^p(r) \left[ \frac{1}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \tilde{\alpha}_p(t,r) \right) \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t)B_i(r) + \tilde{\alpha}_p(t,r) \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t)B'_i(r) \right] \\ + \tilde{\alpha}_p(t,r) \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t)\bar{Z}_i^p(r) \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t)B_i(r) = \left( H \left( \sum_{i=-1}^{N+1} \alpha_i^{(3)}(t)B_i(r) \right) - a_1 \lambda(t) \sum_{i=-1}^{N+1} \alpha_i^{(3)}(t)B_i(r) \right) \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t)B_i(r). \quad (42) $$
For healthy cells $\tilde{S}_h$, similar calculations lead to

$$ \sum_{i=-1}^{N+1} \frac{d\alpha_i^{(2)}}{dt}(t) B_i(r) + \sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) Z_i^s(r) \left[ \frac{1}{r^{d-1}} \frac{\partial}{\partial r} \left( r^{d-1} \tilde{\alpha}_p(t, r) \right) \sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) B_i(r) + \tilde{\alpha}_s(t, r) \sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) B_i'(r) \right] \\ + \tilde{\alpha}_s(t, r) \sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) \bar{Z}_i^s(r) \sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) B_i(r) = -a_2 \lambda(t) \sum_{i=-1}^{N+1} \alpha_i^{(3)}(t) B_i(r) \sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) B_i(r). \quad (43) $$

Similarly, for the third equation describing the concentration of oxygen and nutrients in the tissue, we have

$$ \sum_{i=-1}^{N+1} \frac{d\alpha_i^{(3)}}{dt}(t) B_i(r) - D \left( \frac{d-1}{r} \sum_{i=-1}^{N+1} \alpha_i^{(3)}(t) B_i'(r) + \sum_{i=-1}^{N+1} \alpha_i^{(3)}(t) B_i''(r) \right) = \\ = M_s \sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) B_i(r) \left( 1 - \sum_{i=-1}^{N+1} \alpha_i^{(3)}(t) B_i(r) \right) - \eta \sum_{i=-1}^{N+1} \alpha_i^{(3)}(t) B_i(r) \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t) B_i(r). \quad (44) $$

To determine the values of $\alpha$ at the boundaries, we consider $\tilde{U}_h = (\tilde{P}_h, \tilde{S}_h, \tilde{M}_h)^T$ as natural cubic splines. Natural cubic splines require that their second derivatives vanish at the endpoints of the interval $[0, R]$. Hence, we have:

$$ \frac{\partial^2 \tilde{U}_h}{\partial r^2}(t, r_0) = \begin{pmatrix} \frac{\partial^2 \tilde{P}_h}{\partial r^2}(t, r_0) \\ \frac{\partial^2 \tilde{S}_h}{\partial r^2}(t, r_0) \\ \frac{\partial^2 \tilde{M}_h}{\partial r^2}(t, r_0) \end{pmatrix} = \begin{pmatrix} \frac{1}{h^2} \alpha_{-1}^{(1)}(t) - \frac{2}{h^2} \alpha_0^{(1)}(t) + \frac{1}{h^2} \alpha_1^{(1)}(t) \\ \frac{1}{h^2} \alpha_{-1}^{(2)}(t) - \frac{2}{h^2} \alpha_0^{(2)}(t) + \frac{1}{h^2} \alpha_1^{(2)}(t) \\ \frac{1}{h^2} \alpha_{-1}^{(3)}(t) - \frac{2}{h^2} \alpha_0^{(3)}(t) + \frac{1}{h^2} \alpha_1^{(3)}(t) \end{pmatrix} = \begin{pmatrix} 0, \\ 0, \\ 0. \end{pmatrix} \quad (45) $$

Which leads to

$$ \alpha_{-1}^{j}(t) = 2\alpha_{0}^{j}(t) - \alpha_{1}^{j}(t), \quad \text{for } j = 1, 2, 3. \quad (46) $$

On the other hand, we set the value of the solution in the point $r = r_0$ to be fixed and note it $U_h(t, r_0) := (P_L(t), S_L(t), M_L(t))^T$. Then we have

$$ \tilde{U}_h(t, r_0) = \begin{pmatrix} \tilde{P}_h(t, r_0) \\ \tilde{S}_h(t, r_0) \\ \tilde{M}_h(t, r_0) \end{pmatrix} = \begin{pmatrix} \frac{1}{6} \left( \alpha_{-1}^{(1)}(t) + 4\alpha_0^{(1)}(t) + \alpha_1^{(1)}(t) \right) \\ \frac{1}{6} \left( \alpha_{-1}^{(2)}(t) + 4\alpha_0^{(2)}(t) + \alpha_1^{(2)}(t) \right) \\ \frac{1}{6} \left( \alpha_{-1}^{(3)}(t) + 4\alpha_0^{(3)}(t) + \alpha_1^{(3)}(t) \right) \end{pmatrix} = \begin{pmatrix} \alpha_0^{(1)}(t) \\ \alpha_0^{(2)}(t) \\ \alpha_0^{(3)}(t) \end{pmatrix} = \begin{pmatrix} P_L(t) \\ S_L(t) \\ M_L(t) \end{pmatrix} \quad (47) $$

By an analogue observations, we get for $r = r_{N+1} = R$ the following

$$ \alpha_{N+1}^{j}(t) = 2\alpha_{N}^{j}(t) - \alpha_{N-1}^{j}(t), \quad \text{for } j = 1, 2, 3, \quad (48) $$

and

$$ \begin{pmatrix} \alpha_{N+1}^{(1)}(t) \\ \alpha_{N+1}^{(2)}(t) \\ \alpha_{N+1}^{(3)}(t) \end{pmatrix} = \begin{pmatrix} P_R \\ S_R \\ M_R \end{pmatrix} \quad (49) $$

Let note for $k = p, s$ the following vector valued functions of size $N-1$ as follows:

$$ Z_k(r) := \begin{pmatrix} Z_1^k(r) - Z_{-1}^k(r) \\ Z_2^k(r) \\ \vdots \\ Z_{N-2}^k(r) \\ Z_{N-1}^k(r) - Z_{N+1}^k(r) \end{pmatrix}, \quad \text{and} \quad \bar{Z}_k(r) := \begin{pmatrix} \bar{Z}_1^k(r) - \bar{Z}_{-1}^k(r) \\ \bar{Z}_2^k(r) \\ \vdots \\ \bar{Z}_{N-2}^k(r) \\ \bar{Z}_{N-1}^k(r) - \bar{Z}_{N+1}^k(r) \end{pmatrix}. \quad (50) $$
Furthermore, we define the following functions

$$Q_p(r) := P_L(t) (2Z_{-1}^p(r) + Z_0^p(r)) + P_R (2Z_{N+1}^p(r) + Z_{N+1}^p(r)), \quad (51)$$

$$\overline{Q}_p(r) := P_L(t) (2\overline{Z}_{-1}^p(r) + \overline{Z}_0^p(r)) + P_R (2\overline{Z}_{N+1}^p(r) + \overline{Z}_{N+1}^p(r)), \quad (52)$$

$$Q_s(r) := P_L(t) (2Z_{-1}^s(r) + Z_0^s(r)) + S_R (2Z_{N+1}^s(r) + Z_{N+1}^s(r)), \quad (53)$$

$$\overline{Q}_s(r) := S_L(t) (2\overline{Z}_{-1}^s(r) + \overline{Z}_0^s(r)) + S_R (2\overline{Z}_{N+1}^s(r) + \overline{Z}_{N+1}^s(r)). \quad (54)$$

Then, by using the relations (46), (47), (48) and (49), we can express the nonlocal terms as follows

$$\begin{aligned} \sum_{i=-1}^{N+1} \alpha_i^{(1)}(t) Z_i^p(r) &= \alpha_{-1}^{(1)}(t) Z_{-1}^p(r) + \alpha_0^{(1)}(t) Z_0^p(r) + \sum_{i=1}^{N-1} \alpha_i^{(1)}(t) Z_i^p(r) + \alpha_N^{(1)}(t) Z_N^p(r) + \alpha_{N+1}^{(1)}(t) Z_{N+1}^p(r), \\ &= P_L(t) (2Z_{-1}^p(r) + Z_0^p(r)) + P_R (2Z_{N+1}^p(r) + Z_{N+1}^p(r)) \\ &\quad + \sum_{i=1}^{N-1} \alpha_i^{(1)}(t) Z_i^p(r) - \alpha_1^{(1)}(t) Z_{-1}^p(r) - \alpha_{N-1}^{(1)}(t) Z_{N+1}^p(r), \\ &= Q_p(r) + Z_p^T(r) \alpha^{(1)}(t). \end{aligned}$$

Similarly we have

$$\sum_{i=-1}^{N+1} \alpha_i^{(1)}(t) \overline{Z}_i^p(r) = \overline{Q}_p(r) + \overline{Z}_p^T(r) \alpha^{(1)}(t),$$

$$\sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) Z_i^s(r) = Q_s(r) + Z_s^T(r) \alpha^{(2)}(t),$$

$$\sum_{i=-1}^{N+1} \alpha_i^{(2)}(t) \overline{Z}_i^s(r) = \overline{Q}_s(r) + \overline{Z}_s^T(r) \alpha^{(2)}(t).$$

Now we express $\tilde{U}_h$ at the points $r_i$ for $i = 0, \dots, N$. We start by the equation (42), at $r = r_0$, we obtain

$$\begin{aligned} \frac{d\alpha_0^{(1)}}{dt}(t) = &- (Q_p(r_0) + Z_p^T(r_0)\alpha^{(1)}(t)) \left[ \frac{1}{r_0^{d-1}} \left( \frac{\partial(r^{d-1}\tilde{\alpha}_p)}{\partial r}(t, r_0) \right) P_L(t) + \tilde{\alpha}_p(t, r_0) \frac{1}{h} (\alpha_1^{(1)}(t) - P_L(t)) \right] \\ &- \tilde{\alpha}_p(t, r_0) (\overline{Q}_p(r_0) + \overline{Z}_p^T(r_0)\alpha^{(1)}(t)) P_L(t) + [H(M_L(t)) - a_1\lambda(t)M_L(t)] P_L(t), \end{aligned}$$

$$\begin{aligned} \frac{d\alpha_0^{(2)}}{dt}(t) = &- (Q_s(r_0) + Z_s^T(r_0)\alpha^{(2)}(t)) \left[ \frac{1}{r_0^{d-1}} \left( \frac{\partial(r^{d-1}\tilde{\alpha}_s)}{\partial r}(t, r_0) \right) S_L(t) + \tilde{\alpha}_s(t, r_0) \frac{1}{h} (\alpha_1^{(1)}(t) - S_L(t)) \right] \\ &- \tilde{\alpha}_s(t, r_0) (\overline{Q}_s(r_0) + \overline{Z}_s^T(r_0)\alpha^{(2)}(t)) S_L(t) - a_2\lambda(t)M_L(t)S_L(t), \end{aligned}$$

$$\frac{d\alpha_0^{(3)}}{dt}(t) = \frac{D(d-1)}{rh} (\alpha_1^{(3)}(t) - M_L(t)) + M_s S_L(t) (1 - M_L(t)) - \eta M_L(t) P_L(t). \quad (55)$$

Now, substituting at $r = r_i$, for $i = 1, \dots, N-1$, we get the following

$$\begin{aligned} & \frac{1}{6} \left( \frac{d\alpha_{i-1}^{(1)}}{dt}(t) + 4\frac{d\alpha_i^{(1)}}{dt}(t) + \frac{d\alpha_{i+1}^{(1)}}{dt}(t) \right) = \\ & - (Q_p(r_i) + Z_p^T(r_i)\alpha^{(1)}(t)) \left[ \frac{1}{r_i^{d-1}} \left( \frac{\partial(r^{d-1}\tilde{\alpha}_p)}{\partial r}(t, r_i) \right) \frac{1}{6} (\alpha_{i-1}^{(1)}(t) + 4\alpha_i^{(1)}(t) + \alpha_{i+1}^{(1)}(t)) + \tilde{\alpha}_p(t, r_i) \frac{1}{2h} (\alpha_{i+1}^{(1)}(t) - \alpha_{i-1}^{(1)}(t)) \right] \\ & - \frac{1}{6} \tilde{\alpha}_p(t, r_i) (\overline{Q}_p(r_i) + \overline{Z}_p^T(r_i)\alpha^{(1)}(t)) (\alpha_{i-1}^{(1)}(t) + 4\alpha_i^{(1)}(t) + \alpha_{i+1}^{(1)}(t)) \\ & + \left[ H \left( \frac{1}{6} (\alpha_{i-1}^{(1)}(t) + 4\alpha_i^{(1)}(t) + \alpha_{i+1}^{(1)}(t)) \right) - \frac{1}{6} a_1 \lambda(t) (\alpha_{i-1}^{(3)}(t) + 4\alpha_i^{(3)}(t) + \alpha_{i+1}^{(3)}(t)) \right] \frac{1}{6} (\alpha_{i-1}^{(1)}(t) + 4\alpha_i^{(1)}(t) + \alpha_{i+1}^{(1)}(t)). \end{aligned} \quad (56)$$
$$
\begin{aligned}
\frac{1}{6} \left( \frac{d\alpha_{i-1}^{(2)}}{dt}(t) + 4 \frac{d\alpha_i^{(2)}}{dt}(t) + \frac{d\alpha_{i+1}^{(2)}}{dt}(t) \right) = & - (Q_s(r_i) + Z_s^T(r_i)\alpha^{(2)}(t)) \left[ \frac{1}{r_i^{d-1}} \left( \frac{\partial(r^{d-1}\tilde{\alpha}_s)}{\partial r}(t, r_i) \right) \frac{1}{6} \left( \alpha_{i-1}^{(2)}(t) + 4\alpha_i^{(2)}(t) + \alpha_{i+1}^{(2)}(t) \right) + \tilde{\alpha}_s(t, r_i) \frac{1}{2h} \left( \alpha_{i+1}^{(2)}(t) - \alpha_{i-1}^{(2)}(t) \right) \right] \\
& - \frac{1}{6} \tilde{\alpha}_s(t, r_i) \left( \overline{Q_s}(r_i) + \overline{Z_s}^T(r_i)\alpha^{(2)}(t) \right) \left( \alpha_{i-1}^{(2)}(t) + 4\alpha_i^{(2)}(t) + \alpha_{i+1}^{(2)}(t) \right) \\
& - \frac{1}{36} a_2 \lambda(t) \left( \alpha_{i-1}^{(3)}(t) + 4\alpha_i^{(3)}(t) + \alpha_{i+1}^{(3)}(t) \right) \left( \alpha_{i-1}^{(2)}(t) + 4\alpha_i^{(2)}(t) + \alpha_{i+1}^{(2)}(t) \right).
\end{aligned}
\quad (57)
$$

$$
\begin{aligned}
\frac{1}{6} \left( \frac{d\alpha_{i-1}^{(3)}}{dt}(t) + 4 \frac{d\alpha_i^{(3)}}{dt}(t) + \frac{d\alpha_{i+1}^{(3)}}{dt}(t) \right) = & \\
& D \left( \frac{d-1}{2rh} \left( \alpha_{i+1}^{(3)}(t) - \alpha_{i-1}^{(3)}(t) \right) + \frac{1}{h^2} \left( \alpha_{i-1}^{(3)}(t) - 2\alpha_i^{(3)}(t) + \alpha_{i+1}^{(3)}(t) \right) \right) \\
& + \frac{M_s}{6} \left( \alpha_{i-1}^{(2)}(t) + 4\alpha_i^{(2)}(t) + \alpha_{i+1}^{(2)}(t) \right) \left( 1 - \frac{1}{6} \left( \alpha_{i-1}^{(3)}(t) + 4\alpha_i^{(3)}(t) + \alpha_{i+1}^{(3)}(t) \right) \right) \\
& - \eta \frac{1}{36} \left( \alpha_{i-1}^{(3)}(t) + 4\alpha_i^{(3)}(t) + \alpha_{i+1}^{(3)}(t) \right) \left( \alpha_{i-1}^{(1)}(t) + 4\alpha_i^{(1)}(t) + \alpha_{i+1}^{(1)}(t) \right).
\end{aligned}
\quad (58)
$$

And finally we get for $r = r_N = R$ the following

$$
\begin{align}
\frac{d\alpha_N^{(1)}}{dt}(t) &= - \left( Q_p(r_N) + Z_p^T(r_N)\alpha^{(1)}(t) \right) \left[ \frac{1}{r_N^{d-1}} \left( \frac{\partial(r^{d-1}\tilde{\alpha}_p)}{\partial r}(t, r_N) \right) P_R + \tilde{\alpha}_p(t, r_N) \frac{1}{h} \left( P_R - \alpha_{N-1}^{(1)}(t) \right) \right] \nonumber \\
&\quad - \tilde{\alpha}_p(t, r_N) \left( \overline{Q}_p(r_N) + \overline{Z}_p^T(r_N)\alpha^{(1)}(t) \right) P_R + [H(M_R) - a_1\lambda(t)M_R] P_R, \nonumber \\
\frac{d\alpha_N^{(2)}}{dt}(t) &= - \left( Q_s(r_N) + Z_s^T(r_N)\alpha^{(2)}(t) \right) \left[ \frac{1}{r_N^{d-1}} \left( \frac{\partial(r^{d-1}\tilde{\alpha}_s)}{\partial r}(t, r_N) \right) S_R + \tilde{\alpha}_s(t, r_N) \frac{1}{h} \left( P_R - \alpha_{N-1}^{(1)}(t) \right) \right] \nonumber \\
&\quad - \tilde{\alpha}_s(t, r_N) \left( \overline{Q}_s(r_N) + \overline{Z}_s^T(r_N)\alpha^{(2)}(t) \right) S_R - a_2\lambda(t)M_R S_R, \nonumber \\
\frac{d\alpha_N^{(3)}}{dt}(t) &= \frac{D(d-1)}{rh} \left( M_R - \alpha_{N-1}^{(3)}(t) \right) + M_s S_R(t) (1 - M_R) - \eta M_R P_R. \tag{59}
\end{align}
$$

For the initial time $t = 0$, we have

$$
\frac{1}{6} \begin{pmatrix} 4\alpha_1^{(1)}(0) + \alpha_2^{(1)}(0) \\ \alpha_1^{(1)}(0) + 4\alpha_2^{(1)}(0) + \alpha_3^{(1)}(0) \\ \vdots \\ \alpha_{N-3}^{(1)}(0) + 4\alpha_{N-2}^{(1)}(0) + \alpha_{N-1}^{(1)}(0) \\ \alpha_{N-2}^{(1)}(0) + 4\alpha_{N-1}^{(1)}(0) \\ 4\alpha_1^{(2)}(0) + \alpha_2^{(2)}(0) \\ \alpha_1^{(2)}(0) + 4\alpha_2^{(2)}(0) + \alpha_3^{(2)}(0) \\ \vdots \\ \alpha_{N-3}^{(2)}(0) + 4\alpha_{N-2}^{(2)}(0) + \alpha_{N-1}^{(2)}(0) \\ \alpha_{N-2}^{(2)}(0) + 4\alpha_{N-1}^{(2)}(0) \\ 4\alpha_1^{(3)}(0) + \alpha_2^{(3)}(0) \\ \alpha_1^{(3)}(0) + 4\alpha_2^{(3)}(0) + \alpha_3^{(3)}(0) \\ \vdots \\ \alpha_{N-3}^{(3)}(0) + 4\alpha_{N-2}^{(3)}(0) + \alpha_{N-1}^{(3)}(0) \\ \alpha_{N-2}^{(3)}(0) + 4\alpha_{N-1}^{(3)}(0) \end{pmatrix} = \begin{pmatrix} \tilde{P}_0(r_1) - \frac{1}{6}\tilde{P}_0(r_0) \\ \tilde{P}_0(r_2) \\ \vdots \\ \tilde{P}_0(r_{N-2}) \\ \tilde{P}_0(r_{N-1}) - \frac{1}{6}\tilde{P}_0(r_N) \\ \tilde{S}_0(r_1) - \frac{1}{6}\tilde{S}_0(r_0) \\ \tilde{S}_0(r_2) \\ \vdots \\ \tilde{S}_0(r_{N-2}) \\ \tilde{S}_0(r_{N-1}) - \frac{1}{6}\tilde{S}_0(r_N) \\ \tilde{M}_0(r_1) - \frac{1}{6}\tilde{M}_0(r_0) \\ \tilde{M}_0(r_2) \\ \vdots \\ \tilde{M}_0(r_{N-2}) \\ \tilde{M}_0(r_{N-1}) - \frac{1}{6}\tilde{M}_0(r_N) \end{pmatrix} := \tilde{U}^0. \quad (60)
$$
Which can be expressed as $\mathbb{A}\alpha(0) = \tilde{U}^0$, where $\mathbb{A}$ is a $3(N-1) \times 3(N-1)$ block matrix defined as

$$ \mathbb{A} = \begin{bmatrix} A & 0 & 0 \\ 0 & A & 0 \\ 0 & 0 & A \end{bmatrix}. \qquad (61) $$

Where $\mathbb{A}$ is the following $(N-1) \times (N-1)$ matrix

$$ A = \frac{1}{6} \begin{pmatrix} 4 & 1 & 0 & \cdots & 0 \\ 1 & 4 & 1 & & \vdots \\ 0 & \ddots & \ddots & & \ddots & 0 \\ \vdots & & 1 & 4 & 1 \\ 0 & \cdots & 0 & 1 & 4 \end{pmatrix}. \qquad (62) $$

Furthermore, let define the matrices of size $(N-1) \times (N-1)$ as follows:

$$ A' := \frac{-1}{2h} \begin{pmatrix} 0 & -1 & 0 & \cdots & 0 \\ 1 & 0 & -1 & & \vdots \\ 0 & \ddots & \ddots & & \ddots & 0 \\ \vdots & & 1 & 0 & -1 \\ 0 & \cdots & 0 & 1 & 0 \end{pmatrix}, \quad A'' := \frac{1}{h^2} \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & & \vdots \\ 0 & \ddots & \ddots & & \ddots & 0 \\ \vdots & & 1 & -2 & 1 \\ 0 & \cdots & 0 & 1 & -2 \end{pmatrix}, \\ \mathbb{M}_Z^k := [\mathbb{Z}_k(r_1), \dots, \mathbb{Z}_k(r_{N-1})]^T, \quad \text{and} \quad \bar{\mathbb{M}}_Z^k := [\bar{\mathbb{Z}}_k(r_1), \dots, \bar{\mathbb{Z}}_k(r_{N-1})]^T, \quad k = p, s. \qquad (63) $$

Let also define for $k = p, s$, the following vectors of size $N-1$ as follows:

$$ \begin{align*} \mathbb{M}_Q^k &:= (Q_k(r_1), \dots, Q_k(r_{N-1}))^T, & \bar{\mathbb{M}}_Q^k &:= (\bar{Q}_k(r_1), \dots, \bar{Q}_k(r_{N-1}))^T, \\ \mathbb{M}_{\tilde{\alpha}}^k(t) &:= (\tilde{\alpha}_k(t, r_1), \dots, \tilde{\alpha}_k(t, r_{N-1}))^T, & \bar{\mathbb{M}}_{\tilde{\alpha}}^k(t) &:= \left( \frac{1}{r_1^{d-1}} \frac{\partial}{\partial r} (r^{d-1} \tilde{\alpha}_p(t, r_1)), \dots, \frac{1}{r_{N-1}^{d-1}} \frac{\partial}{\partial r} (r^{d-1} \tilde{\alpha}_p(t, r_{N-1})) \right)^T, \\ V_M(t) &= \left( \frac{M_L(t)}{6}, 0, \dots, 0, \frac{M_R}{6} \right)^T, & \bar{r}_d &:= (d-1) \left( \frac{1}{r_1^{d-1}}, \dots, \frac{1}{r_{N-1}^{d-1}} \right)^T, \quad \text{and} \quad 1 = (1, \dots, 1)^T. \tag{64} \end{align*} $$

Subsequently, $\alpha$ satisfies the following ordinary differential equation (ODE):

$$ \left\{ \begin{aligned} \mathbb{A} \frac{d\alpha}{dt}(t) &= \Phi(t, \alpha(t)), & t \in [0, T], \\ \mathbb{A} \alpha(0) &= \tilde{U}^0. \end{aligned} \right. \qquad (65) $$

Where $\Phi$ is a nonlinear function defined by the equation (55),(56),(57),(58) and (59). Then, by using the notations (61),(63) and (64), $\Phi$ may be written as follows:

$$ \Phi(t, \alpha(t)) = - \begin{pmatrix} \mathbb{M}_Z^p & 0 & 0 \\ 0 & \mathbb{M}_Z^s & 0 \\ 0 & 0 & 0 \end{pmatrix} \alpha(t) + \begin{bmatrix} \mathbb{M}_Q^p \\ \mathbb{M}_Q^s \\ 0 \end{bmatrix} \circ \begin{pmatrix} \overline{\mathbb{M}}_{\tilde{\alpha}}^p(t) \\ \overline{\mathbb{M}}_{\tilde{\alpha}}^s(t) \\ 0 \end{pmatrix} \circ \mathbb{A} + \begin{bmatrix} \mathbb{M}_{\tilde{\alpha}}^p \\ \mathbb{M}_{\tilde{\alpha}}^s \\ 0 \end{bmatrix} \circ \begin{bmatrix} A' & 0 & 0 \\ 0 & A' & 0 \\ 0 & 0 & 0 \end{bmatrix} \alpha(t) \\ + \begin{bmatrix} (H(A\alpha^{(3)}(t) + V_M(t)) - a_1\lambda(t)A\alpha^{(3)}(t)) \circ (A\alpha^{(1)}(t)) \\ -a_2\lambda(t)(A\alpha^{(3)}(t)) \circ (A\alpha^{(1)}(t)) \\ D(\bar{r}_d \odot A' + A'')\alpha^{(3)}(t) + M_s(A\alpha^{(2)}(t)) \circ (1 - A\alpha^{(3)}(t)) - \eta(A\alpha^{(3)}(t)) \circ (A\alpha^{(1)}(t)) \end{bmatrix} + f(t, \alpha(t)). $$

Where $\odot$ denotes for the Hadamard product, and $f$ is a $3 \times (N-1)$ vector, expressed as follows:

$$ \begin{aligned} f_1(t, \alpha(t)) = \frac{-1}{6} \frac{d\alpha_0^{(1)}}{dt}(t) - (Q_p(r_1) + \mathbb{Z}_p^T(r_1)\alpha^{(1)}(t)) & \left[ \frac{1}{6r_1^{d-1}} \left( \frac{\partial(r^{d-1}\tilde{\alpha}_p)}{\partial r}(t, r_1) \right) P_L(t) - \tilde{\alpha}_p(t, r_1) \frac{1}{2h} P_L(t) \right] \\ & - \frac{1}{6} \tilde{\alpha}_p(t, r_1) \left( \overline{Q}_p(r_1) + \overline{\mathbb{Z}}_p^T(r_1)\alpha^{(1)}(t) \right) P_L(t) \\ & - \frac{1}{36} a_1 \lambda(t) \left[ M_L(t) (P_L(t) + 4\alpha_1^{(1)}(t) + \alpha_2^{(1)}(t)) + P_L(t) (M_L(t) + 4\alpha_1^{(3)}(t) + \alpha_2^{(3)}(t)) \right], \end{aligned} $$
$$
\begin{aligned}
f_{N-1}(t, \alpha(t)) = \frac{-1}{6} \frac{d\alpha_N^{(1)}}{dt}(t) - (Q_p(r_{N-1}) + Z_p^T(r_{N-1})\alpha^{(1)}(t)) & \left[ \frac{1}{6r_{N-1}^{d-1}} \left( \frac{\partial(r^{d-1}\widetilde{\alpha}_p)}{\partial r}(t, r_{N-1}) \right) P_R + \widetilde{\alpha}_p(t, r_{N-1}) \frac{1}{2h} P_R \right] \\
& - \frac{1}{6}\widetilde{\alpha}_p(t, r_{N-1}) (\overline{Q}_p(r_{N-1}) + \overline{Z}_p^T(r_{N-1})\alpha^{(1)}(t)) P_R \\
& - \frac{1}{36} a_1 \lambda(t) \left[ M_R \left( P_R + 4\alpha_{N-1}^{(1)}(t) + \alpha_{N-2}^{(1)}(t) \right) + P_R \left( M_R + 4\alpha_{N-1}^{(3)}(t) + \alpha_{N-2}^{(3)}(t) \right) \right],
\end{aligned}
$$

$$
\begin{aligned}
f_N(t, \alpha(t)) = \frac{-1}{6} \frac{d\alpha_0^{(2)}}{dt}(t) - (Q_s(r_1) + Z_s^T(r_1)\alpha^{(2)}(t)) & \left[ \frac{1}{6r_1^{d-1}} \left( \frac{\partial(r^{d-1}\widetilde{\alpha}_s)}{\partial r}(t, r_1) \right) S_L(t) - \widetilde{\alpha}_s(t, r_1) \frac{1}{2h} S_L(t) \right] \\
& - \frac{1}{6}\widetilde{\alpha}_s(t, r_1) (\overline{Q}_s(r_1) + \overline{Z}_s^T(r_1)\alpha^{(2)}(t)) S_L(t) \\
& - \frac{1}{36} a_2 \lambda(t) \left[ M_L(t) (S_L(t) + 4\alpha_1^{(2)}(t) + \alpha_2^{(2)}(t)) + S_L(t) (M_L(t) + 4\alpha_1^{(3)}(t) + \alpha_2^{(3)}(t)) \right],
\end{aligned}
$$

$$
\begin{aligned}
f_{2N-2}(t, \alpha(t)) = \frac{-1}{6} \frac{d\alpha_N^{(2)}}{dt}(t) - (Q_S(r_{N-1}) + Z_S^T(r_{N-1})\alpha^{(2)}(t)) & \left[ \frac{1}{6r_{N-1}^{d-1}} \left( \frac{\partial(r^{d-1}\widetilde{\alpha}_S)}{\partial r}(t, r_{N-1}) \right) S_R + \widetilde{\alpha}_S(t, r_{N-1}) \frac{1}{2h} S_R \right] \\
& - \frac{1}{6}\widetilde{\alpha}_S(t, r_{N-1}) (\overline{Q}_S(r_{N-1}) + \overline{Z}_S^T(r_{N-1})\alpha^{(2)}(t)) S_R \\
& - \frac{1}{36} a_2 \lambda(t) \left[ M_R \left( S_R + 4\alpha_{N-1}^{(2)}(t) + \alpha_{N-2}^{(2)}(t) \right) + S_R \left( M_R + 4\alpha_{N-1}^{(3)}(t) + \alpha_{N-2}^{(3)}(t) \right) \right],
\end{aligned}
$$

$$
\begin{aligned}
f_{2N-1}(t, \alpha(t)) = \frac{-1}{6} \frac{d\alpha_0^{(3)}}{dt}(t) + D \left( \frac{1}{h^2} M_L(t) - \frac{d-1}{2rh} M_L(t) \right) + \frac{M_s}{6} S_L(t) - \frac{M_s}{36} \left[ M_L(t) \left( S_L(t) + 4\alpha_1^{(2)}(t) + \alpha_2^{(2)}(t) \right) \right. \\
\left. + S_L(t) \left( M_L(t) + 4\alpha_1^{(3)}(t) + \alpha_2^{(3)}(t) \right) \right] - \frac{\eta}{36} \left[ M_L(t) \left( P_L(t) + 4\alpha_1^{(1)}(t) + \alpha_2^{(1)}(t) \right) + P_L(t) \left( M_L(t) + 4\alpha_1^{(3)}(t) + \alpha_2^{(3)}(t) \right) \right],
\end{aligned}
$$

$$
\begin{aligned}
f_{3N-3}(t, \alpha(t)) = \frac{-1}{6} \frac{d\alpha_0^{(3)}}{dt}(t) + D \left( \frac{1}{h^2} M_R + \frac{d-1}{2rh} M_R \right) + \frac{M_s}{6} S_R - \frac{M_s}{36} \left[ M_R \left( S_R + 4\alpha_{N-1}^{(2)}(t) + \alpha_{N-2}^{(2)}(t) \right) \right. \\
\left. + S_R \left( M_L(t) + 4\alpha_{N-1}^{(3)}(t) + \alpha_{N-2}^{(3)}(t) \right) \right] - \frac{\eta}{36} \left[ M_R \left( P_L(t) + 4\alpha_{N-1}^{(1)}(t) + \alpha_{N-2}^{(1)}(t) \right) + P_R \left( M_R + 4\alpha_{N-1}^{(3)}(t) + \alpha_{N-2}^{(3)}(t) \right) \right],
\end{aligned}
$$

and $f_j(t, \alpha(t)) = 0$ for all $j \notin \{1, N-1, N, 2N-2, 2N-1, 3N-3\}$.

The solution of the ODE (65) can be approached using several methods found in the literature. However, due to the stiffness of the equation, which is a common characteristic in systems derived from discretized partial differential equations.

Common explicit methods, such as the popular Runge-Kutta family of integrators, are not suitable for such stiff problems. To overcome this limitation, we employ the Backward Differentiation Formula (BDF) method, which is an implicit multi-step method specifically designed for solving stiff ordinary differential equations. This method allows for larger time steps while maintaining accuracy and stability, which can be done by using information from previous time steps and solving an implicit equation at each step. Here, we apply the BDF method of order $p = 6$, for maximal stability and efficiency. The method applied to (65) gives the following approximations

$$ A \alpha_n - \beta \Delta t \Phi(t_n, \alpha_n) - A \sum_{j=0}^{p-1} \eta_j \alpha_{n-j-1} = 0. \qquad (66) $$

In this context, $\Delta t = \frac{T}{p}$ represents the time step, and $\alpha_n = [\alpha_1^n, \dots, \alpha_{N-1}^n]^T$ is the BDF method's approximation of the vector $\alpha(t)$ at time $t_n = n\Delta t$ for $n = 0, \dots, p$. The coefficients $\eta_j$ and $\beta$ for the $p$-step BDF formula are provided in Table 2.
<table><caption>Table 2: BDF coefficients for order 6.</caption><thead><tr><th>p</th><th>β</th><th>η<sub>0</sub></th><th>η<sub>1</sub></th><th>η<sub>2</sub></th><th>η<sub>3</sub></th><th>η<sub>4</sub></th><th>η<sub>5</sub></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td><sup>2</sup>&frasl;<sub>3</sub></td><td><sup>4</sup>&frasl;<sub>3</sub></td><td><sup>-1</sup>&frasl;<sub>3</sub></td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td><sup>6</sup>&frasl;<sub>11</sub></td><td><sup>18</sup>&frasl;<sub>11</sub></td><td><sup>-9</sup>&frasl;<sub>11</sub></td><td><sup>2</sup>&frasl;<sub>11</sub></td><td></td><td></td><td></td></tr><tr><td>4</td><td><sup>12</sup>&frasl;<sub>25</sub></td><td><sup>48</sup>&frasl;<sub>25</sub></td><td><sup>-36</sup>&frasl;<sub>25</sub></td><td><sup>16</sup>&frasl;<sub>25</sub></td><td><sup>-3</sup>&frasl;<sub>25</sub></td><td></td><td></td></tr><tr><td>5</td><td><sup>60</sup>&frasl;<sub>137</sub></td><td><sup>300</sup>&frasl;<sub>137</sub></td><td><sup>-300</sup>&frasl;<sub>137</sub></td><td><sup>200</sup>&frasl;<sub>137</sub></td><td><sup>-75</sup>&frasl;<sub>137</sub></td><td><sup>12</sup>&frasl;<sub>137</sub></td><td></td></tr><tr><td>6</td><td><sup>60</sup>&frasl;<sub>147</sub></td><td><sup>360</sup>&frasl;<sub>147</sub></td><td><sup>-450</sup>&frasl;<sub>147</sub></td><td><sup>400</sup>&frasl;<sub>147</sub></td><td><sup>-225</sup>&frasl;<sub>147</sub></td><td><sup>72</sup>&frasl;<sub>147</sub></td><td><sup>-10</sup>&frasl;<sub>147</sub></td></tr></tbody></table>

At each time step $n$, the equation (66) for $\alpha_n$ must be solved by reformulating it into the following form:

$$ \varphi(\alpha_n) = A\alpha_n - \beta\Delta t\Phi(t_n, \alpha_n) - A\sum_{j=0}^{p-1} \eta_j \alpha_{n-j-1} = 0. \quad (67) $$

To solve the equation (66) for $\alpha_n$ efficiently, we employ the Newton method. We used the solution from the previous time step as an initial guess. The Newton method approximates $\alpha_n$ through a series of iterations ($\zeta_k$) as follows:

$$ \begin{cases} \zeta_0, \\ \zeta_{k+1} = \zeta_k - [J_{\varphi}(\zeta_k)]^{-1} \varphi(\zeta_k), \quad k \ge 0, \end{cases} \quad (68) $$

where $J_{\varphi}(\zeta_k)$ is the Jacobian matrix of $\varphi$, expressed as follows:

$$ J_{\varphi}(\zeta_k) = A - \beta \Delta t J_{\Phi}(\zeta_k), \quad (69) $$

with $J_{\Phi}$ denotes for the Jacobian matrix of $\Phi$ with respect to $\alpha_n$.

## 5. Numerical experiments

In this section, we aim to examine the accuracy of our proposed numerical scheme by comparing it with a known analytical solution. This evaluation is essential to determine the efficiency and validity of the method. The construction of a classical solution for a simplified version of the system presented in equation (34) presents significant challenges, mainly due to the nonlocal term, which must be calculated analytically in order to prevent errors that could arise from numerical integration. Then in section 5.2, we use values similar to those reported by Lefebvre in [4], examining the dynamics of avascular tumor growth in the absence of medical treatment, considering a nonlocal velocity.

Throughout this section, we consider the spatial dimension of our system (10) to be $d = 2$. In this context, the nonlocal terms are expressed as follows:

$$ Z_i^k(r) := (\tilde{\gamma}_p * B_i)(r) = 2 \int_0^R \int_0^\pi \tilde{\gamma}_p \left( \sqrt{s^2 + r^2 + 2rs \cos(\theta)} \right) B(s) \, ds \, d\theta \quad \text{and} $$

$$ \bar{Z}_i^k(r) := \frac{\partial (\tilde{\gamma}_p * B_i)(r)}{\partial r} = 2 \int_0^R \int_0^\pi \frac{r + s \cos(\theta)}{\sqrt{r^2 + s^2 + 2rs \cos(\theta)}} \tilde{\gamma}_p' \left( \sqrt{s^2 + r^2 + 2rs \cos(\theta)} \right) B(s) \, ds \, d\theta. $$

These double integrals are calculated numerically using the function integral2 in MATLAB, which is effective for this type of two-dimensional integration over rectangular regions.

The simulations performed in section 5.1 are executed in the spatial domain $[0, 1]$, while in section 5.2, the spatial domain $[10^{-3}, 2]$ is considered. It is important to note that radial hyperbolic and radial parabolic equations are likely to suffer from singularities at the origin (see [32]). To mitigate this issue, the calculations are initiated from a small radius in space. The numerical simulations presented in this paper were performed using MATLAB on a standard laptop computer equipped with an Intel Core i5 processor (2.6 GHz, 2 cores) and 8GB of RAM. No specialized GPU acceleration was utilized for these computations. This configuration was sufficient for running our simulations, demonstrating the efficiency of the proposed numerical method.
## 5.1. Convergence test with known analytical solution

In this section, we consider the following coupled system with source terms

$$ 
\left\{
\begin{array}{l}
\displaystyle \frac{\partial \tilde{P}}{\partial t}(t,r) + 2\pi r \frac{\partial}{\partial r} \left( \left( \int_0^R \tilde{P}(t,s)ds \right) \tilde{P} \right)(t,r) = \left( H(\tilde{M})(t,r) - t\tilde{M}(t,r) - 2 \left( \int_0^R \tilde{P}(t,s)ds \right) \right) \tilde{P}(t,r) + f_1(t,r), \\[1.5em]
\displaystyle \frac{\partial \tilde{S}}{\partial t}(t,r) + 2\pi r \frac{\partial}{\partial r} \left( \left( \int_0^R \tilde{S}(t,s)ds \right) \tilde{S} \right)(t,r) = - \left( t\tilde{M}(t,r) + 2 \left( \int_0^R \tilde{S}(t,s)ds \right) \right) \tilde{S}(t,r) + f_2(t,r), \\[1.5em]
\displaystyle \frac{\partial \tilde{M}}{\partial t}(t,r) - 0.01 \left( r \frac{\partial^2 \tilde{M}}{\partial r^2}(t,r) + \frac{\partial \tilde{M}}{\partial r}(t,r) \right) = \tilde{S}(t,x)(1 - \tilde{M}(t,x)) - \tilde{M}(t,x)\tilde{P}(t,x) + f_3(t,r), \\[1.5em]
\tilde{U}(t,R) = \tilde{U}_{ex}(t,R), \quad \tilde{U}(t,0) = \tilde{U}_{ex}(t,0), \\[1.5em]
\tilde{U}(0,r) = \tilde{U}_{ex}(0,r),
\end{array}
\right.
\quad (t,r) \in [0, T] \times [0, R], \qquad (70) 
$$

The source term $(f(t,r) = f_1(t,r), f_2(t,r), f_3(t,r))$ is chosen in order that the system (70) has the following analytical solution

$$ \tilde{U}_{ex}(t,r) = \begin{pmatrix} \tilde{P}_{ex}(t,r) \\ \tilde{S}_{ex}(t,r) \\ \tilde{M}_{ex}(t,r) \end{pmatrix} = \begin{pmatrix} \frac{3}{2} \exp\left(\frac{\left(r - \frac{t}{2000}\right)^2}{50}\right) + \exp\left(\frac{\left(r - \frac{t}{350}\right)^2}{100}\right) + \frac{1}{2} \exp\left(\frac{\left(r - \frac{t}{400}\right)^2}{150}\right) \\ \frac{3}{2} \exp\left(\frac{\left(r - \frac{t}{2000}\right)^2}{70}\right) + \exp\left(\frac{\left(r - \frac{t}{350}\right)^2}{140}\right) + \frac{1}{2} \exp\left(\frac{\left(r - \frac{t}{400}\right)^2}{210}\right) \\ \frac{3}{2} \exp\left(\frac{\left(r - \frac{t}{2000}\right)^2}{80}\right) + \frac{4}{5} \exp\left(\frac{\left(r - \frac{t}{200}\right)^2}{160}\right) + \frac{3}{5} \exp\left(\frac{\left(r - \frac{t}{300}\right)^2}{240}\right) \end{pmatrix}, \quad (71) $$

The nonlocal system (70) can be considered a particular case of the general system (34) by setting specific parameter values. Specifically, we choose $\tilde{\alpha}_p(t,r) = \tilde{\alpha}_s(t,r) = r$, $\tilde{\gamma}_p(r) = \tilde{\gamma}_s(r) = 1$, $a_1 = 1$, $a_2 = 1$, $\lambda(t) = t$, $M_s = \eta = 1$, and $D = \frac{r}{100}$. These parameters were selected to allow the radial convolution to be calculated analytically. This is important because radial systems often exhibit singularities at $r = 0$. To address these singularities, the velocity functions $\tilde{\alpha}_p$, $\tilde{\alpha}_s$, and the diffusion parameter $D$ are specially chosen.

To numerically solve the system, we discretize the domain $[0, R]$ using a uniform grid. The grid points are defined as $r_i = ih$, for $i = 1, \dots, N$. Where $h$ is the spatial step size given by $h = \frac{R}{N}$, and $N$ being the number of discretization nodes.

We then compare the approximate numerical solution to the analytical solution using the relative error metrics defined as follows:

$$ E_n^\infty = \frac{\left\| \tilde{U}(t^n) - \tilde{U}_h(t^n) \right\|_\infty}{\left\| \tilde{U}(t^n) \right\|_\infty} = \frac{\sup_{0 \le i \le N} \left\| \tilde{U}(t^n, r_i) - \tilde{U}_h(t^n, r_i) \right\|}{\sup_{0 \le i \le N} \left\| \tilde{U}(t^n, r_i) \right\|}, $$

$$ \text{and} \quad E^\infty = \frac{\left\| \tilde{U} - \tilde{U}_h \right\|_\infty}{\left\| \tilde{U} \right\|_\infty} = \frac{\sup_n \left\| \tilde{U}(t^n) - \tilde{U}_h(t^n) \right\|_\infty}{\sup_n \left\| \tilde{U}(t^n) \right\|_\infty}. $$

These relative error metrics provide a means to quantitatively assess the accuracy of the numerical solution compared to the analytical solution. The following table presents the execution times and relative errors for different numbers of nodes:

<table><caption>**Table 3:** Execution time and relative error for different node numbers for $R = 1$, and $T = 10$.</caption><tr><td>N</td><td>50</td><td>100</td><td>150</td><td>200</td><td>250</td><td>300</td></tr><tr><td>E<sup>∞</sup></td><td>3.344e-03</td><td>1.456e-03</td><td>6.610e-04</td><td>3.707e-04</td><td>2.440e-04</td><td>1.649e-04</td></tr><tr><td>CPU Time</td><td>18.96s</td><td>48.56s</td><td>188.53s</td><td>422.28s</td><td>853.93s</td><td>1619.54s</td></tr></table>
The following figure illustrates the comparison between the approximate solution $\tilde{U}_h$ and the exact solution $\tilde{U}_{ex}$ of our nonlocal system at various time points: $t = 0, t = 25, t = 50, t = 75$, and $t = 100$.

<figure><img src="image_4.png" /><figcaption>Figure 1: Approximate solution $\tilde{U}_h$ and Exact solution $\tilde{U}_{ex}$ for $t = 0, t = 25, t = 50, t = 75, t = 100$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Approximation of $\tilde{M}$ on $[0, 1] \times [0, 100]$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2: Approximation of $\tilde{P}$ on $[0, 1] \times [0, 100]$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3: Approximation of $\tilde{S}$ on $[0, 1] \times [0, 100]$.</figcaption></figure>
For a range of time values, a comparison between the numerical solution, denoted as $\tilde{U}_a$, and the exact solution, denoted as $\tilde{U}_{ex}$, is conducted, as illustrated in Figure 1. In this figure, it is evident that the numerical solution approximates the exact solution with high accuracy across the different time instances (t=0, t=25, t=50, t=75, t=100). This demonstrates the effectiveness of the numerical method in closely replicating the exact solution over the specified time intervals.

<table><caption>**Table 4:** Relative error $E^\infty$.</caption><thead><tr><th>N</th><th>100</th><th>200</th><th>300</th><th>400</th><th>500</th></tr></thead><tbody><tr><th>h</th><td>0.01</td><td>0.005</td><td>0.0033</td><td>0.0025</td><td>0.002</td></tr><tr><th>Δt</th><td>0.005</td><td>0.0025</td><td>0.00165</td><td>0.00125</td><td>0.001</td></tr><tr><th>t=20</th><td>6.9801e-04</td><td>1.7519e-04</td><td>7.7218e-05</td><td>4.3550e-05</td><td>2.7949e-05</td></tr><tr><th>t=40</th><td>4.1530e-04</td><td>1.0432e-04</td><td>5.3297e-05</td><td>2.6102e-05</td><td>1.67076e-05</td></tr><tr><th>t=60</th><td>8.0758e-04</td><td>2.0209e-04</td><td>9.0763e-05</td><td>4.9907e-05</td><td>3.2346e-05</td></tr><tr><th>t=80</th><td>1.4171e-03</td><td>3.5577e-04</td><td>1.5894e-04</td><td>8.9162e-05</td><td>5.6989e-05</td></tr><tr><th>t=100</th><td>1.8271e-03</td><td>4.5819e-04</td><td>2.0402e-04</td><td>1.1462e-04</td><td>7.3364e-05</td></tr><tr><th>E<sup>∞</sup></th><td>5.0599e-04</td><td>1.2579e-04</td><td>7.0742e-05</td><td>3.1597e-05</td><td>2.0176e-05</td></tr></tbody></table>

Table 4 presents various values of the relative error $E_n^\infty$ for $\Delta t = \frac{h}{2}$ with different spatial step sizes $h$. The data suggest that as $h$ decreases, the computational error also decreases. Although proving the theoretical convergence of the Backward Differentiation Formula (BDF) method for a nonlinear Ordinary Differential Equation (ODE) is challenging, the computed results indicate a reduction in error with smaller spatial steps. This observation is consistent with existing literature on numerical analysis and the stability of numerical integration methods. The known stability of the BDF method in stiff problems corroborates the observed error reduction trends with decreasing $h$ values.

## 5.2. Numerical Test for Tumor Growth

In this subsection, we present numerical simulations of the nonlocal tumor growth model (9), under the scenario of no medical intervention.

The initial conditions for our simulations are defined as follows:

$$ \tilde{P}(0, r) = \frac{0.3e^2}{e^2 + e^{20r}}, \quad \tilde{S}(0, r) = \frac{0.8e^2 + e^{20r}}{e^2 + e^{20r}}, \quad \tilde{M}(0, r) = 1. $$

At the boundary radius $R$, the conditions are specified as follows:

$$ \tilde{P}(t, R) = 0, \quad \tilde{S}(t, R) = 1, \quad \tilde{M}(t, R) = 1. $$

The velocities and kernel functions incorporated in our model are expressed as follows.

$$ \tilde{\alpha}_p(t, r) = \frac{2e}{e + e^{10r}}, \quad \tilde{\alpha}_s(t, r) = \frac{0.2e^{7.5}}{e^{7.5} + e^{15r}}, \quad \tilde{\gamma}_p(r) = 0.2 \exp(-\frac{r}{50}), \quad \tilde{\gamma}_s(r) = \frac{1}{100} \exp(-\frac{r}{150}). $$

We used a hypoxia threshold value of $M_{th} = 0.5$. The growth factor parameter, $\kappa$, is set at 0.02. The model also integrates a diffusion coefficient $D$ of 0.8. Furthermore, the nutrient consumption rate, $\eta$, is 1.5, and the blood vessel creation rate is characterized by $M_s = 1$.

The results of our numerical simulations are demonstrated in Figure 5.
<figure><img src="image_8.png" /><figcaption>(a) $t = 5$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) $t = 20$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) $t = 40$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(d) $t = 60$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(e) $t = 80$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(f) $t = 100$</figcaption></figure>

**Figure 5:** Radial tumor cells density $\tilde{P}$, healthy cells density $\tilde{S}$, and the concentration of nutriment and oxygen $\tilde{M}$ for different time instances.

<figure><img src="image_14.png" /><figcaption>(a) $t = 0$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) $t = 20$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) $t = 40$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) $t = 60$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(e) $t = 80$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(f) $t = 100$</figcaption></figure>

**Figure 6:** Tumor cells density $P$ for different time instances.

Figures 6 and 7 illustrates the growth of the tumor core, particularly at its edges where higher nutrient concentrations are present, facilitating cellular proliferation. The presence of abundant nutrients at the periphery enables
<figure><img src="image_20.png" /><figcaption>Figure 7: Radial tumor cells density for different time instances.</figcaption></figure>

the cells to rapidly divide and expand outward. In contrast, the center of the tumor exhibits a decrease in tumor cell density due to insufficient oxygen supply, leading to necrosis. The lack of oxygen, known as hypoxia, results in cell death and the formation of a necrotic core. This process highlights the dynamic nature of tumor growth, where peripheral areas continue to expand while central regions undergo cell death. Furthermore, the hypoxic conditions in the tumor center can lead to the release of signaling molecules that may promote angiogenesis or new blood vessel formation in an attempt to supply the tumor with more oxygen and nutrients.

## Conclusion

The process of tumor growth remains a complex subject to study, predominantly due to the multifaceted nature of the factors that influence it. A detailed understanding of the mechanisms governing tumor growth is essential for the advancement of effective anticancer strategies.

In this paper, we have presented an innovative continuum model of tumor growth, which includes both cell-cell and cell-matrix adhesion and interactions. This model is characterized by the integration of nonlocal terms within a system of partial differential equations (PDEs). We have successfully demonstrated the existence and uniqueness of solutions through the application of semigroup properties and Banach's fixed-point theorem, in combination with results established by Keimer in [22]. Furthermore, we have developed a numerical scheme for resolving the radial nonlocal system. This scheme employs B-splines for spatial discretization and Backward Differentiation Formula (BDF) for temporal discretization. Extensive numerical tests have been conducted to validate the model's accuracy, and we have demonstrated experimental results related to tumor growth.

The use of nonlocal balance equations in our model offers several significant advantages over traditional local approaches. For example, nonlocal models can avoid the formation of shock waves and singularities that often challenge local hyperbolic systems. This regularizing effect leads to smoother solutions and may prevent the breakdown of numerical schemes, which is a common issue in local models of complex phenomena. The efficient numerical solution of nonlocal systems has potential applications in many fields such as population dynamics, material science, and climate modeling, where long-range effects play a crucial role and where avoiding discontinuities is critical.

In the future, our research aims to optimize the factors of medical drugs within the system (9), using the theory of optimal control for nonlocal balance equations. This will help design new strategies that minimize the auxiliary effects of treatment and maximize its efficacy. Using mathematical models, we intend to precisely adjust drug dosages, timings, and combinations to achieve the best therapeutic outcomes while reducing negative side effects.

## Acknowledgment

This work was supported by the Centre National pour la Recherche Scientifique et Technique (CNRST) of Morocco and the Laboratoire de Mathématiques Pures et Appliquées (LMPA) at the Université du Littoral Côte d'Opale (ULCO), France.
Bibliography

[1] Aziz Belmiloudi. Mathematical modeling and optimal control problems in brain tumor targeted drug delivery strategies. *International Journal of Biomathematics*, 10(04):1750056, May 2017.

[2] Luigi Preziosi, editor. *Cancer modelling and simulation*. Chapman & Hall/CRC mathematical biology and medicine series. Chapman & Hall/CRC, Boca Raton, Fla, 2003.

[3] Helen Byrne and Luigi Preziosi. Modelling solid tumour growth using the theory of mixtures. *Mathematical medicine and biology: a journal of the IMA*, 20(4):341–366, 2003.

[4] Guillaume Lefebvre, François Cornelis, Patricio Cumsille, Thierry Colin, Clair Poignard, and Olivier Saut. Spatial modelling of tumour drug resistance: the case of gist liver metastases. *Mathematical medicine and biology: a journal of the IMA*, 34(2):151–176, 2017.

[5] Didier Bresch, Thierry Colin, Emmanuel Grenier, Benjamin Ribba, and Olivier Saut. Computational Modeling of Solid Tumor Growth: The Avascular Stage. *SIAM Journal on Scientific Computing*, 32(4):2321–2344, January 2010.

[6] Amina Eladdadi, Peter Kim, and Dann Mallet, editors. *Mathematical Models of Tumor-Immune System Dynamics*, volume 107 of *Springer Proceedings in Mathematics & Statistics*. Springer New York, New York, NY, 2014.

[7] Tamal Das, Kai Safferling, Sebastian Rausch, Niels Grabe, Heike Boehm, and Joachim P. Spatz. A molecular mechanotransduction pathway regulates collective migration of epithelial cells. *Nature Cell Biology*, 17(3):276–287, March 2015.

[8] Thomas S Deisboeck and Iain D Couzin. Collective behavior in cancer cell populations. *Bioessays*, 31(2):190–197, 2009.

[9] C.T Lee, M.F Hoopes, J Diehl, W Gilliland, G Huxel, E.V Leaver, K Mccann, J Umbanhowar, and A Mogilner. Non-local Concepts and Models in Biology. *Journal of Theoretical Biology*, 210(2):201–219, May 2001.

[10] Dieter Armbruster, Pierre Degond, and Christian Ringhofer. A model for the dynamics of large queuing networks and supply chains. *SIAM Journal on Applied Mathematics*, 66(3):896–920, 2006.

[11] Rinaldo M Colombo and Elena Rossi. Hyperbolic predators vs. parabolic prey (vol 13, pg 369, 2015). *COMMUNICATIONS IN MATHEMATICAL SCIENCES*, 13(5):1345–1346, 2015.

[12] Rinaldo M. Colombo, Mauro Garavello, and Magali Lécureux-Mercier. A CLASS OF NONLOCAL MODELS FOR PEDESTRIAN TRAFFIC. *Mathematical Models and Methods in Applied Sciences*, 22(04):1150023, April 2012.

[13] Aekta Aggarwal and Paola Goatin. Crowd dynamics through non-local conservation laws. *Bulletin of the Brazilian Mathematical Society, New Series*, 47(1):37–50, March 2016.

[14] Dieter Armbruster, Daniel E Marthaler, Christian Ringhofer, Karl Kempf, and Tae-Chang Jo. A continuum model for a re-entrant factory. *Operations research*, 54(5):933–950, 2006.

[15] Hao Kang, Shigui Ruan, and Xiao Yu. Age-Structured Population Dynamics with Nonlocal Diffusion. *Journal of Dynamics and Differential Equations*, 34(2):789–823, June 2022.

[16] Ariel Ramírez-Torres, Salvatore Di Stefano, and Alfio Grillo. Influence of non-local diffusion in avascular tumour growth. *Mathematics and Mechanics of Solids*, 26(9):1264–1293, 2021.

[17] Marvin Fritz, Ernesto A. B. F. Lima, Vanja Nikolić, J. Tinsley Oden, and Barbara Wohlmuth. Local and nonlocal phase-field models of tumor growth and invasion due to ECM degradation. *Mathematical Models and Methods in Applied Sciences*, 29(13):2433–2468, December 2019. arXiv:1906.07788 [math].
[18] Marvin Fritz, Ernesto A. B. F. Lima, J. Tinsley Oden, and Barbara Wohlmuth. On the unsteady Darcy–Forchheimer–Brinkman equation in local and nonlocal tumor growth models. *Mathematical Models and Methods in Applied Sciences*, 29(09):1691–1731, August 2019.

[19] A. Gerisch and M.A.J. Chaplain. Mathematical modelling of cancer cell invasion of tissue: Local and non-local models and the effect of adhesion. *Journal of Theoretical Biology*, 250(4):684–704, February 2008.

[20] Vasiliki Bitsouni and Raluca Eftimie. Non-local Parabolic and Hyperbolic Models for Cell Polarisation in Heterogeneous Cancer Cell Populations. *Bulletin of Mathematical Biology*, 80(10):2600–2632, October 2018.

[21] Alexander Keimer and Lukas Pflug. Existence, uniqueness and regularity results on nonlocal balance laws. *Journal of Differential Equations*, 263(7):4023–4069, October 2017.

[22] Alexander Keimer, Lukas Pflug, and Michele Spinola. Existence, uniqueness and regularity of multi-dimensional nonlocal balance laws with damping. *Journal of Mathematical Analysis and Applications*, 466(1):18–55, October 2018.

[23] Helen M. Byrne. A weakly nonlinear analysis of a model of avascular solid tumour growth. *Journal of Mathematical Biology*, 39(1):59–89, July 1999.

[24] H.M. Byrne, J.R. King, D.L.S. McElwain, and L. Preziosi. A two-phase model of solid tumour growth. *Applied Mathematics Letters*, 16(4):567–573, May 2003.

[25] Thomas Michel, Jérôme Fehrenbach, Valérie Lobjois, Jennifer Laurent, Aurelie Gomes, Thierry Colin, and Clair Poignard. Mathematical modeling of the proliferation gradient in multicellular tumor spheroids. *Journal of theoretical biology*, 458:133–147, 2018.

[26] Vittorio Cristini, Xiangrong Li, John S Lowengrub, and Steven M Wise. Nonlinear simulations of solid tumor growth using a mixture model: invasion and branching. *Journal of mathematical biology*, 58:723–763, 2009.

[27] T. Michel, J. Fehrenbach, V. Lobjois, J. Laurent, A. Gomes, T. Colin, and C. Poignard. Mathematical modeling of the proliferation gradient in multicellular tumor spheroids. *Journal of Theoretical Biology*, 458:133–147, December 2018.

[28] Songmu Zheng. *Nonlinear evolution equations*. Chapman and Hall/CRC, 2004.

[29] Thierry Cazenave and Alain Haraux. *An introduction to semilinear evolution equations*, volume 13. Oxford University Press, 1998.

[30] Haim Brezis. *Functional Analysis, Sobolev Spaces and Partial Differential Equations*. Universitext. Springer, New York, 2010.

[31] Alexander Keimer, Lukas Pflug, and Michele Spinola. Nonlocal scalar conservation laws on bounded domains and applications in traffic flow. *SIAM Journal on Mathematical Analysis*, 50(6):6271–6306, 2018.

[32] Kamran Mohseni and Tim Colonius. Numerical treatment of polar coordinate singularities. *Journal of Computational Physics*, 157(2):787–795, 2000.