# Mathematical investigation of normal and abnormal wound healing dynamics: local and non-local models

O E Adebayo, S. Urcun, G. Rolin, S P A Bordas, D. Trucu, Raluca Eftimie

► To cite this version:

O E Adebayo, S. Urcun, G. Rolin, S P A Bordas, D. Trucu, et al., Mathematical investigation of normal and abnormal wound healing dynamics: local and non-local models. Mathematical Biosciences and Engineering, 2023, 20 (9), pp.17446-17498. 10.3934/mbe.2023776. hal-04353465

HAL Id: hal-04353465

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

Submitted on 19 Dec 2023

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Research article

# Mathematical investigation of normal and abnormal wound healing dynamics: local and non-local models

O. E. Adebayo<fnref n="1" />, S. Urcun<fnref n="2" />, G. Rolin<fnref n="3" />,<fnref n="4" />, S. P. A. Bordas<fnref n="2" />, D. Trucu<fnref n="5" />, R. Eftimie<fnref n="1" />,<fnref n="5" />,*

**Abstract:** The movement of cells during (normal and abnormal) wound healing is the result of bio-mechanical interactions that combine cell responses to growth factors as well as cell-cell and cell-matrix interactions (adhesion and remodelling). It is known that cells can communicate and interact locally and non-locally with other cells inside the tissues, through mechanical forces that act locally and at a distance, as well as through long non-conventional cell protrusions. In this study, we consider a non-local PDE model for the interactions between fibroblasts, macrophages, and extracellular matrix via a growth factor ($TGF-\beta$) in the context of wound healing. For the non-local interactions, we consider two types of kernels (a Gaussian kernel and a cone-shaped kernel), two types of cell-ECM adhesion functions (adhesion only to higher ECM densities vs. adhesion to higher/lower densities) and two types of cell proliferation terms (with and without decay due to overcrowding). We investigate numerically the dynamics of this non-local model, as well as the dynamics of the localised versions of this model (obtained when cell perception radius reduces towards 0). The results suggest that: (i) local models explain normal wound healing and non-local models could explain also abnormal wound healing (although the results are parameter-dependent); (ii) the models can explain two types of wound healing (by primary intention, when the wounds margins come together from the side; and by secondary intention when the wound healing from bottom up).

**Keywords:** normal and abnormal wound healing; non-local models; local models; FEM;

<footnote n="1">¹ Laboratoire de mathématiques de Besançon, UMR CNRS 6623, Université de Franche-Comté, Besançon 25000, France</footnote>

<footnote n="2">² Department of Engineering, Faculty of Science, Technology and Medicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg</footnote>

<footnote n="3">³ INSERM CIC-1431, CHU Besançon, Besançon 25000, France</footnote>

<footnote n="4">⁴ INSERM, EFS BFC, UMR1098, RIGHT Interactions, Greffon-Hôte-Tumeur/Ingénierie Cellulaire et Génique, Université Franche-Comté 25000, France</footnote>

<footnote n="5">⁵ Division of Mathematics, University of Dundee, Dundee, DD1 4HN, United Kingdom</footnote>

<footnote>* Correspondence: E-mail: raluca.eftimie@univ-fcomte.fr, r.a.eftimie@dundee.ac.uk.</footnote>
# 1. Introduction

Wound healing is a complex process through which an organism tries to restore the integrity of biological tissues following their physical damage [1,2]. While different tissues in the body are able to heal after wounding, here we focus only on normal and abnormal wounds associated with skin tissue. There are two main ways for cutaneous tissue to heal: (i) by primary intention, where the incision is narrow and wound heals as its edges are brought close together; (ii) by secondary intention, where the wound heals from the bottom of the wound up. A normal wound healing (either by primary or secondary intention) follows a series of interconnected phases (hemostasis, inflammation, proliferation, and remodelling) [3], which eventually returns the wounded tissue close to its original state. However, when there are dysregulation during the various phases of wound healing, it can lead to abnormal wounds and excessive scars [4], such as hypertrophic scars and keloids. In regard to these fibrotic pathologies, we first note that both hypertrophic and keloid scars can occur following healing by primary and/or secondary intention. Second, we note that although both hypertrophic scars and keloids are both characterised by the presence of excessive scar tissue, only keloids grow beyond the borders of the original wound in a tumour-like manner [5,6].

Unfortunately, neither hypertrophic scars nor keloids benefit from satisfactory prevention tools or treatment, and thus remain an uncovered clinical need. Deciphering the overlapping phases involved in normal and abnormal wound healing is key to understanding how these fibrosis are triggered and evolve overtime. In the following, we briefly detail key biological aspects associated with each of these four phases:

*   **Hemostasis:** following an injury, blood fills the wound area and blood platelets coagulate to form a fibrin mesh (i.e., blood clot) that prevents further blood loss. Moreover, this fibrin mesh (which forms the basis of the new extracellular matrix inside the wounded tissue) acts as a scaffold for early cell migration into the wound. Platelets also release different inflammatory cytokines and growth factors, such as TGF-$\beta$, which promote the inflammatory phase [7–9]. In fact, as we will see below, TGF-$\beta$ has a critical role in the different phases of wound healing [10], and for this reason throughout this study we focus on this growth factor.

*   **Inflammation:** this phase is characterised by increased capillary permeability and cell infiltration and migration into the wound site [7]. The neutrophils that arrive within hours to the wound area to clean the debris, are followed within 3-4 days by macrophages that also deterse the wound and also secrete various cytokines and growth factors, such as TGF-$\beta$, which then recruit fibroblasts and initiate the formation of granulation tissue [7,9].

*   **Proliferation:** this phase is characterised by ECM proteins deposition and connective tissue contraction through the differentiation of fibroblasts into myofibroblasts, and finally wound re-epithelization. During proliferation phase, fibroblasts are chemo-attracted and migrate (namely in response to TGF-$\beta$), proliferate and produce ECM components (e.g., collagens and fibronectin) that lead to the formation of granulation tissue [7,9].

*   **Remodelling:** this phase is characterised by a re-organisation of collagen fibres in granulation tissue [9], with type III collagen being replaced with newly secreted type-I collagen. Such remodelling is the result of a "synthesis-degradation" balance under the control of matrix metalloproteinases (MMPs) that cause collagen breakdown. These proteinases, and their inhibitors (TIMPs) are secreted by fibroblasts and macrophages as well. TGF-$\beta$ can inhibit the secretion
of these metalloproteinases leading to accumulation of collagen fibres [7] as observed in both hypertrophic scars and keloids.

Since the complexity of the above wound healing processes makes it challenging to interpret the experimental results, over the past three decades various mathematical models have been derived to investigate the biological mechanisms involved in wound healing [11, 12]. Some of these models focus on normal wounds [13–16], while others focus on abnormal wounds [17–19]. Models usually describe healing by first intention [20], focus mainly on bio-chemical interactions between cells (or between cells and ECM) and do not include bio-mechanical forces involved in wound contraction. In contrast, models that focus on mechanical interactions [2, 9, 21, 22] usually describe healing by second intention [20]. However, to our knowledge, there are not many models that investigate at the same time the biological mechanisms behind healing by first or secondary intentions, in the context of understanding normal vs. abnormal healing phenomena. Also, there are not many models that investigate keloid formation as a result of abnormal wound healing (see the review in [12]).

The goal of this study is to shed new light on the biological mechanisms that could generate different types of wounds (normal or abnormal), that heal by primary or by secondary intention. To this end, we consider a theoretical (i.e., modelling and numerical) approach to investigate some of the inflammatory and bio-mechanical aspects involved in the different phases of wound healing discussed above. Here we focus on the roles of TGF-β, macrophages, and fibroblasts during the inflammation, proliferation, and remodelling steps. The non-local model developed here accounts for the non-local impact of bio-mechanical interactions between cells, and among cells and extracellular matrix (ECM)(as recent experimental studies have shown long-distance interactions between cells [23–25]). By considering different non-local interaction kernels, different cell adhesion functions and cell proliferation functions, we aim to shed light on the potential mechanisms involved in normal and abnormal wounds healed by first or second intent.

The paper is structured as follows. In Section 2, we describe the non-local model and present its reduction to a local model when we assume that the interactions become very localised. In Section 3, we describe briefly the finite element approach used to discretise the equations and present a series of numerical simulations of the local and non-local models. We conclude in Section 4 with a summary and discussion of the results

## 2. Description of the mathematical model

To describe the wound healing process, we focus on the coupled dynamics of the following main variables: the concentration of a growth factor such as TGF-β, g(x, t), the density of fibroblasts, f(x, t), the density of macrophages, m(x, t), and the density of the ECM, e(x, t). For simplicity, throughout this study, we use the compact vector notation $u = (g, f, m, e)^T$. The time and space evolution of these variables is described by the following equations:

$$ \frac{\partial g}{\partial t} = D_g \Delta g - \lambda_g g + W(f, m), \qquad (2.1a) $$

$$ \frac{\partial f}{\partial t} = \nabla \cdot (D_f \nabla f - \mu_f f A_f [g, f, m, e]) - \lambda_f f + p_f(g)f(1 - \rho(u)), \qquad (2.1b) $$

$$ \frac{\partial m}{\partial t} = \nabla \cdot (D_m \nabla m - \mu_m m A_m [g, f, m, e]) - \lambda_m m + p_m(g)m(1 - \rho(u)), \qquad (2.1c) $$
$$ \frac{\partial e}{\partial t} = -e(\alpha_f f + \alpha_m m) + p_e(f) e(1 - \rho(u)). \qquad (2.1d) $$

These equations incorporate the following biological assumptions:

* The changes in the concentration of the growth factor (Eq (2.1a)) at a position $(\mathbf{x}, t)$, are the result of the diffusion of these molecules (with diffusion coefficient $D_g$) [26], the decay of the growth factor (at a rate $\lambda_g$) [7, 27] and the secretion of this growth factor by fibroblasts and macrophages [28–30], as described by the term $W(f, m)$ below:

$$ W(f, m) = p_g^f f + p_g^m m, \qquad (2.2) $$

where $p_g^f$ is the production rate of the growth factor by fibroblasts [28–30], and $p_g^m$ is the production rate of the growth factor by macrophages [28].

* The changes in the density of fibroblasts (Eq (2.1b)) at a position $(\mathbf{x}, t)$, are the result of the flux of the fibroblasts [28, 31] (which consists of a linear diffusion term with coefficient $D_f$, and a cell migration term with coefficient $\mu_f$, which is a consequence of cell-cell and cell-ECM adhesion [32, 33]). In addition, the change in fibroblasts densities is the result of fibroblasts’ apoptosis [28] (at a rate $\lambda_f$), and fibroblasts self-renewal via proliferation [31, 34] (at a rate $p_f(g)$ which depends on the concentration of the growth factor $g$). Note that Since various experimental studies have shown that cells can interact at a distance with other cells [23–25], in Eq (2.1b), $A_f[g, f, m, e]$ denotes a non-local spatial flux term describing the adhesion processes between the fibroblasts, macrophages, and ECM, responsible for the directed movement of the fibroblasts (and the role of growth factor in these processes). A detailed description of this non-local flux term is given below. Moreover,

$$ \rho(u) = w_g g + w_f f + w_m m + w_e e, $$

is the cumulative volume fraction space occupied by the components of our system. Here $w_g, w_f, w_m, w_e > 0$ are indices for the volume fraction spaces occupied by the growth factor, fibroblasts, macrophages, and the ECM respectively. Therefore, the proliferation of cells is described by the logistic term $f(1 - \rho(u))$ that models cell proliferation under nutrient availability, which is consistent with some experimental studies [35].

* The changes in the density of macrophages (Eq (2.1c)) at a position $(\mathbf{x}, t)$, are a result of the flux of macrophages [36] (which consists of a linear diffusion term with coefficient $D_m$, a cell migration term with coefficient $\mu_m$, which is a consequence of cell-cell and cell-ECM adhesion [32, 33]), macrophages’ apoptosis (at a rate $\lambda_m$), and the logistic proliferation of macrophages (at a rate $p_m(g)$ which depends on the growth factor $g$) [37]. Here $A_m[g, f, m, e]$ denotes a non-local spatial flux term describing the adhesion processes between the fibroblasts, macrophages, ECM (and the role of TGF-$\beta$ on these adhesive interactions). The proliferation of macrophages is described by the logistic term $m(1 - \rho(u))$.

* The changes in the density of ECM (Eq (2.1d)) are the result of degradation [38, 39] caused by ECM degrading enzymes secreted by fibroblasts [40] at a constant rate $\alpha_f$, and by macrophages [41, 42] at a constant rate $\alpha_m$. It is also the result of ECM remodelling [37, 43] at a rate $p_e(f)$ assumed to depend linearly on fibroblasts’ density (since dermal fibroblasts are closely linked to ECM as both producer and resident cells [44]). Finally, as in [45, 46], we consider a logistic-type remodelling of ECM.
The adhesive cell-cell and cell-ECM interactions between the cells distributed at **x** and the surrounding cells and ECM perceived over a ball-shaped sensing region **B**(**x**, **R**) := **x** + **B**(0, **R**) of radius **R** > 0 are expressed via the non-local terms:

$$A_{f,m}[g, f, m, e](\mathbf{x}, t) = \frac{1}{R} \int_{\mathbf{B}(0,R)} K(\|\mathbf{y}\|_2) \mathbf{n}(\mathbf{y})(1 - \rho(u))^{+} \Gamma_{f,m}(\mathbf{x} + \mathbf{y}, t) \, dy. \quad (2.3)$$

Here, **B**(0, **R**) := {$\zeta \in \mathbb{R}^2 : \|\zeta\|_2 \le R$} is the usual closed ball of radius **R** centred at 0, and **n**(**y**) denotes the unit radial vector originating from **x** and moving towards **x** + **y** ∈ **B**(0, **R**) for any **y** ∈ **B**(0, **R**) and is defined as follows:

$$\mathbf{n}(\mathbf{y}) := \begin{cases} \frac{\mathbf{y}}{\|\mathbf{y}\|_2}, & \text{if } \mathbf{y} \in \mathbf{B}(0, R) \setminus \{(0, 0)\} \\ (0, 0), & \text{otherwise,} \end{cases} \quad (2.4)$$

where $\|\cdot\|_2$ is the usual Euclidean norm. The kernel $K(\cdot) : [0, R] \to [0, 1]$ is a radially symmetric kernel that gives the interaction range of cells (i.e., interactions between the reference cell at position **x** and neighbours at **x** + **y**). Examples of such kernels are

a. Gaussian kernel (see Figure 1(a))

$$K_1(z) = \frac{1}{2\pi\sigma^2} e^{-\frac{z^2}{2\sigma^2}}. \quad (2.5)$$

b. Cone-shaped kernel (see Figure 1(b))

$$K_2(z) = \frac{3}{\pi R^2} \left(1 - \frac{z}{R}\right). \quad (2.6)$$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Kernels describing the long-distance cell-cell and cell-ECM interactions. (a) Gaussian kernel (see Eq (2.5)) with standard deviation $\sigma = 0.04$; (b) cone-shaped kernel (see Eq (2.6)) with $R = 0.1$. In Eqs (2.6) and (2.5), $z = \sqrt{y_1^2 + y_2^2}$.</figcaption></figure>

Moreover, the term $(1 - \rho(u))^+ := \max\{(1 - \rho(u)), 0\}$ is included to avoid overcrowding within the non-local interactions. Finally, the function $\Gamma_{f,m}(\mathbf{x} + \mathbf{y}, t)$ describes the type and magnitude of cell-cell and
cell-ECM adhesive interactions between cells at position **x** and neighbours at position **x** + **y**. To define these functions $\Gamma_{f,m}$, we assume that fibroblasts are co-cultured with macrophages on ECM. Note here that some experimental studies observed that macrophages cannot adhere to ECM (at least to type-I collagen) [47], while other experimental studies showed that macrophages can adhere to some type of substrate (e.g., cross-linked hydrogel) [48]. To investigate these contradictory experimental results, in this study, we consider two sub-cases:

I. No macrophage-ECM adhesion.

$$ \Gamma_f(\mathbf{x} + \mathbf{y}, t) := \mathbf{S}_{ff} f(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{fm} m(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{fe} e(\mathbf{x} + \mathbf{y}, t), \quad (2.7a) $$

$$ \Gamma_m(\mathbf{x} + \mathbf{y}, t) := \mathbf{S}_{mm} m(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{mf} f(\mathbf{x} + \mathbf{y}, t). \quad (2.7b) $$

II. Including macrophage-ECM adhesion ($S_{me} > 0$ in $\Gamma_m$).

$$ \Gamma_f(\mathbf{x} + \mathbf{y}, t) := \mathbf{S}_{ff} f(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{fm} m(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{fe} e(\mathbf{x} + \mathbf{y}, t), \quad (2.8a) $$

$$ \Gamma_m(\mathbf{x} + \mathbf{y}, t) := \mathbf{S}_{mm} m(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{mf} f(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{me} e(\mathbf{x} + \mathbf{y}, t). \quad (2.8b) $$

In the above equations, we considered the strengths of fibroblast-fibroblast interactions ($S_{ff}$), fibroblast-macrophage interactions ($S_{fm}$), fibroblast-ECM interactions ($S_{fe}$), macrophage-macrophage interactions ($S_{mm}$), macrophage-fibroblast interactions ($S_{mf}$), macrophage-ECM interactions ($S_{me}$). Since these interaction strengths depend on TGF-$\beta$ and on the presence of ECM [49], to define them we consider a monotonically-increasing Hill-type function (depending on “$e + g$” and satisfying $S_j(e, g)=0$ for $e=0$ and $g=0$):

$$ S_j := S_j^{\max} \frac{e + g}{1 + e + g}, \quad j \in \{ff, fm, mm, mf, fe, me\}, \quad (2.9) $$

The above function is used to describe cell movement up-gradient of ECM and neighbouring cells in the presence of the growth factor (and will be used throughout this study). However, in the early stages of wound healing, the fibroblasts and macrophages leave the healthy tissue to move onto the newly-formed fibrin mesh at the bottom of the wound, to help with the formation of granulation tissue. To describe this dynamics, in Section 3.3, we use the following function for cell-cell and cell-ECM adhesion:

$$ S_j := S_j^{\max} \frac{(e - e_c) + g}{1 + e + g}, \quad j \in \{ff, fm, mm, mf, fe, me\}, \quad (2.10) $$

where $e_c$ is a matrix threshold for the transition between up-gradient cell movement and down-gradient cell movement.

Because we do not really know the spatial range over which this adhesive cell-cell and cell-ECM interactions have an impact (although some experimental studies have suggested that cells can sense up to a few rows of neighbouring cell [50]), in the next section we also consider localised versions of model (2.1). To this end, we assume that the cell perception radius $R \to 0$ and thus we can use classical Taylor expansions to transform the non-local interactions into local interactions.
## 2.1. Reduction of the non-local models to local models: Gaussian kernel

The non-local term in Eq (2.3), with the Gaussian kernel defined by Eq (2.5) on a square domain (see also Figure 1(a)) can be written as:

$$A_{(\circ)}[g, f, m, e](\mathbf{x}, t) := \frac{1}{R} \int_{-R}^{R} \int_{-R}^{R} \mathbf{K}(\mathbf{y})\mathbf{n}(\mathbf{y})(G(u(\mathbf{x} + \mathbf{y}, t))) dy_1 dy_2, \quad (2.11)$$

where $A_{(\circ)}[g, f, m, e] \in \{A_m[g, f, m, e], A_f[g, f, m, e]\}$, $\mathbf{y} = (y_1, y_2)$. The function $G(u(\mathbf{x} + \mathbf{y}, t))$ is defined by the following expressions:

* for $A_f$ we have

$$G(u(\mathbf{x} + \mathbf{y}, t)) = (\mathbf{S}_{ff}f(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{fm}m(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{fe}e(\mathbf{x} + \mathbf{y}, t))(1 - \rho(u))^+ \quad (2.12)$$

* for $A_m$ we have the following general expression

$$G(u(\mathbf{x} + \mathbf{y}, t)) = (\mathbf{S}_{mm}m(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{mf}f(\mathbf{x} + \mathbf{y}, t) + \mathbf{S}_{me}e(\mathbf{x} + \mathbf{y}, t))(1 - \rho(u))^+ \quad (2.13)$$

Case I. (i.e., Eq (2.7), no macrophage-ECM adhesion) corresponds to $S_{me} = 0$ in Eq (2.13), while Case II. (i.e., Eq (2.8), with macrophage-ECM adhesion) corresponds to $S_{me} > 0$.

The truncated local Taylor expansion of $G(u(\mathbf{x} + \mathbf{y}, t))$ around $\mathbf{x} := (x_1, x_2)$ is given as

$$G(u(\mathbf{x} + \mathbf{y}, t)) = G(u(\mathbf{x}, t)) + y_1 \frac{\partial}{\partial x_1} G(u(\mathbf{x}, t)) + y_2 \frac{\partial}{\partial x_2} G(u(\mathbf{x}, t)) + O(y_1^2 + y_2^2). \quad (2.14)$$

Inserting Eq (2.14) into Eq (2.11) we have that

$$A_{(\circ)}[g, f, m, e](\mathbf{x}, t) = \frac{1}{R} \int_{-R}^{R} \int_{-R}^{R} \mathbf{K}(\mathbf{y})\mathbf{n}(\mathbf{y}) \left( G(u(\mathbf{x}, t)) + y_1 \frac{\partial}{\partial x_1} G(u(\mathbf{x}, t)) + y_2 \frac{\partial}{\partial x_2} G(u(\mathbf{x}, t)) + O(y_1^2 + y_2^2) \right) dy_1 dy_2. \quad (2.15)$$

Choosing a kernel as in Eq (2.5), and $\mathbf{n}(\mathbf{y})$ defined in Eq (2.4), and inserting them into Eq (2.15) (ignoring the higher order terms) we have that:

$$A_{(\circ)}[g, f, m, e](\mathbf{x}, t) = \frac{1}{R} G(u(\mathbf{x}, t)) \left( \int_{-R}^{R} \int_{-R}^{R} \frac{y_1}{2\pi\sigma^2\sqrt{y_1^2 + y_1^2}} \exp\left(-\frac{y_1^2 + y_2^2}{2\sigma^2}\right) dy_1 dy_2, \right. \\ \left. \int_{-R}^{R} \int_{-R}^{R} \frac{y_2}{2\pi\sigma^2\sqrt{y_1^2 + y_1^2}} \exp\left(-\frac{y_1^2 + y_2^2}{2\sigma^2}\right) dy_1 dy_2 \right) +$$
$$ 
\begin{aligned}
&+\frac{1}{R}\left(\frac{\partial}{\partial x_1}G(u(\mathbf{x}, t))\right) \int_{-R}^{R} \int_{-R}^{R} \frac{y_1^2}{2\pi\sigma^2\sqrt{y_1^2 + y_1^2}} \exp\left(-\frac{y_1^2 + y_2^2}{2\sigma^2}\right) dy_1 dy_2, \\
&\frac{\partial}{\partial x_2}G(u(\mathbf{x}, t)) \int_{-R}^{R} \int_{-R}^{R} \frac{y_2^2}{2\pi\sigma^2\sqrt{y_1^2 + y_1^2}} \exp\left(-\frac{y_1^2 + y_2^2}{2\sigma^2}\right) dy_1 dy_2.
\end{aligned} 
\quad (2.16) $$

Following the steps detailed in Appendix 4, we obtain that $A_{(\circ)}[g, f, m, e] \to 0$ as $R \to 0$ for all considered cases (see Eq (4.1)). Hence, model (2.1) reduces to:

$$ \frac{\partial g}{\partial t} = D_g \Delta g - \lambda_g g + W(f, m), \quad (2.17a) $$

$$ \frac{\partial f}{\partial t} = \nabla \cdot (D_f \nabla f) - \lambda_f f + p_f(g)f(1 - \rho(u)), \quad (2.17b) $$

$$ \frac{\partial m}{\partial t} = \nabla \cdot (D_m \nabla m) - \lambda_m m + p_m(g)m(1 - \rho(u)), \quad (2.17c) $$

$$ \frac{\partial e}{\partial t} = -e(\alpha_f f + \alpha_m m) + p_e(f) e(1 - \rho(u)). \quad (2.17d) $$

We emphasize that the difference between the localised model (2.17) and the non-local model (2.1) (for the case of the Gaussian kernel), is that cell motility in the local model is given only by the diffusion terms.

## 2.2. Reduction of the non-local models to local models: cone-shaped kernel.

Now we consider as in Gerisch and Chaplain [51], a cone-shaped kernel given by Eq (2.6) on a circular domain with a radius $R$. Then, the non-local term in Eq (2.3) takes the form:

$$ A_{(\circ)}[g, f, m, e](\mathbf{x}, t) := \frac{1}{R} \int_{0}^{R} r \int_{0}^{2\pi} \mathbf{n}(\theta) K(r) G(u(\mathbf{x} + r\mathbf{n}(\theta), t)) d\theta dr. \quad (2.18) $$

Here, $\mathbf{n}(\theta) = (\cos \theta, \sin \theta)^T$ is a vector denoting the outer unit normal of the angle $\theta$.

* For the fibroblasts non-local flux term $A_f[g, f, m, e]$, the cell-cell and cell-ECM interactions are given by

$$ G(u(\mathbf{x} + r\mathbf{n}(\theta), t)) = (\mathbf{S}_{ff} f(\mathbf{x} + r\mathbf{n}(\theta), t) + \mathbf{S}_{fm} m(\mathbf{x} + r\mathbf{n}(\theta), t) + \mathbf{S}_{fe} e(\mathbf{x} + r\mathbf{n}(\theta), t))(1 - \rho(u)). $$

* For the macrophages non-local flux term $A_m[g, f, m, e]$, the cell-cell and cell-ECM interactions are given by

$$ G(u(\mathbf{x} + r\mathbf{n}(\theta), t)) = (\mathbf{S}_{mm} m(\mathbf{x} + r\mathbf{n}(\theta), t) + \mathbf{S}_{mf} f(\mathbf{x} + r\mathbf{n}(\theta), t) + \mathbf{S}_{me} e(\mathbf{x} + r\mathbf{n}(\theta), t))(1 - \rho(u)). $$

As before, in the above equation condition $\mathbf{S}_{me} = 0$ corresponds to Case I. (no macrophage-ECM adhesion), while condition $\mathbf{S}_{me} > 0$ corresponds to Case II. (macrophage-ECM adhesion present).
The Taylor series expansion of $G(u(\mathbf{x} + r\mathbf{n}(\theta), t)) = (G \circ u)(\mathbf{x} + r\mathbf{n}(\theta), t)$ around $r = 0$ is given as:

$$G(u(\mathbf{x} + r\mathbf{n}(\theta), t)) = G(u(\mathbf{x}, t)) + \left\langle \frac{d}{dr} G(u(\mathbf{x}, t)), r\mathbf{n}(\theta) \right\rangle + \dots \qquad (2.19)$$

where

$$\frac{d}{dr} G(u(\mathbf{x}, t)) = \nabla_u G(u(\mathbf{x}, t)) \nabla u(\mathbf{x}, t).$$

Assuming that the functions $G$ and $u$ are smooth, Eq (2.18) becomes

$$A_{(\circ)}[g, f, m, e] = \frac{1}{R} \int_{0}^{R} r \int_{0}^{2\pi} \mathbf{n}(\theta) K(r) G(u(\mathbf{x}, t)) d\theta dr \\ + \frac{1}{R} \int_{0}^{R} r \int_{0}^{2\pi} \mathbf{n}(\theta) K(r) \langle \nabla_u G(u(\mathbf{x}, t)) \nabla u(\mathbf{x}, t), r\mathbf{n}(\theta) \rangle d\theta dr. \qquad (2.20)$$

For $K(r)$ chosen as in Eq (2.6) and using the expression in Eq (4.2) in Appendix (4), we obtain that for $R \to 0$, $A_{(\circ)}[g, f, m, e] \to A_{(\circ)}^{0}(g, f, m, e)$. Since the limit function $A_{(\circ)}^{0}$ depends on $G$, in the following, we describe in detail this limit function. To this end, we also assume that $(1 - \rho(u)) > 0$. Moreover, we write $u(\mathbf{x}, t)$ as $u$ to simplify the notation.

In the limit $R \to 0$, the non-local flux terms $A_{fm}[g, f, m, e]$ approach the local terms $A_{fm}^{0}(g, f, m, e)$, where

$$A_{f}^{0}(g, f, m, e) = \frac{1}{4} \left\{ (1 - \rho(u)) \left( \mathbf{S}_{ff} \nabla f + \mathbf{S}_{fm} \nabla m + \mathbf{S}_{fe} \nabla e + f \left( \frac{\partial}{\partial g} \mathbf{S}_{ff} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{ff} \nabla e \right) \right. \right. \\ + m \left( \frac{\partial}{\partial g} \mathbf{S}_{fm} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{fm} \nabla e \right) + e \left( \frac{\partial}{\partial g} \mathbf{S}_{fe} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{fe} \nabla e \right) \\ \left. - \left( \mathbf{S}_{ff} f + \mathbf{S}_{fm} m + \mathbf{S}_{fe} e \right) \left( w_{g} \nabla g + w_{f} \nabla f + w_{m} \nabla m + w_{e} \nabla e \right) \right\}, \qquad (2.21a)$$

$$A_{m}^{0}(g, f, m, e) = \frac{1}{4} \left\{ (1 - \rho(u)) \left( \mathbf{S}_{mf} \nabla f + \mathbf{S}_{mm} \nabla m + \mathbf{S}_{me} \nabla e + f \left( \frac{\partial}{\partial g} \mathbf{S}_{mf} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{mf} \nabla e \right) \right) \right. \\ + m \left( \frac{\partial}{\partial g} \mathbf{S}_{mm} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{mm} \nabla e \right) + e \left( \frac{\partial}{\partial g} \mathbf{S}_{me} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{me} \nabla e \right) \\ \left. - \left( \mathbf{S}_{mf} f + \mathbf{S}_{mm} m + \mathbf{S}_{fe} e \right) \left( w_{g} \nabla g + w_{f} \nabla f + w_{m} \nabla m + w_{e} \nabla e \right) \right\}. \qquad (2.21b)$$

In Eq (2.21a), the term $\mathbf{S}_{ff} \nabla f + \mathbf{S}_{fm} \nabla m + \mathbf{S}_{fe} \nabla e$ on the right-hand side describes fluxes due to fibroblast-fibroblast, fibroblast-macrophage, and fibroblast-ECM adhesions, with strengths $\mathbf{S}_{ff}$, $\mathbf{S}_{fm}$, and $\mathbf{S}_{fe}$ respectively. The term $f(\frac{\partial}{\partial g} \mathbf{S}_{ff} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{ff} \nabla e) + m(\frac{\partial}{\partial g} \mathbf{S}_{fm} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{fm} \nabla e) + e(\frac{\partial}{\partial g} \mathbf{S}_{fe} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{fe} \nabla e)$ describes the change in the strength of fibroblast-fibroblast, fibroblast-macrophage, and fibroblast-ECM adhesion with respect to the growth factor and ECM, along the positive gradients of growth factor and ECM respectively. Finally, the term $-(\mathbf{S}_{ff} f + \mathbf{S}_{fm} m + \mathbf{S}_{fe} e)(w_g \nabla g + w_f \nabla f + w_m \nabla m + w_e \nabla e)$
prevents the uncontrolled accumulation of fibroblasts, macrophages, and ECM by directing them down the gradients of growth factor, fibroblasts, macrophages and ECM.

In Eq (2.21b), the term $S_{mf} \nabla f + S_{mm} \nabla m$ on the right-hand side represents macrophage-fibroblast and macrophage-macrophage adhesion with strengths $S_{mf}$ and $S_{mm}$ respectively. The term $f(\frac{\partial}{\partial g}S_{mf} \nabla g + \frac{\partial}{\partial e}S_{mf} \nabla e) + m(\frac{\partial}{\partial g}S_{mm} \nabla g + \frac{\partial}{\partial e}S_{mm} \nabla e)$ represents the change in the macrophage-fibroblast and macrophage-macrophage adhesion strength with respect to growth factor and ECM, up the positive gradients of growth factor and ECM, respectively. The term $-(S_{mf}f + S_{mm}m)(w_g \nabla g + w_f \nabla f + w_m \nabla m + w_e \nabla e)$ prevents the uncontrolled accumulation of macrophages and fibroblasts by directing them down the negative gradients of growth factors, fibroblasts, macrophages and ECM.

The difference between Case I. (no macrophage-ECM adhesion; $S_{me} = 0$) and Case II. (macrophage-ECM adhesion present; $S_{me} > 0$) consists in the additional terms on the right-hand side of Eq (2.21b) that contain $S_{me}$. The term $e(\frac{\partial}{\partial g}S_{me} \nabla g + \frac{\partial}{\partial e}S_{me} \nabla e)$ represents the change in the macrophage-ECM adhesion strength with respect to the growth-factor and ECM, along the positive gradients of growth-factor and ECM respectively.

The term $-S_{me}e(w_g \nabla g + w_f \nabla f + w_m \nabla m + w_e \nabla e)$ prevents the uncontrolled accumulation of macrophages-ECM, by directing them down the negative gradients of the growth-factor, fibroblasts, macrophages, and the ECM.

Substituting Eq (2.21) into Eq (2.1) leads to the following local model:

$$ \frac{\partial g}{\partial t} = D_g \Delta g - \lambda_g g + W(f, m), \qquad (2.22a) $$

$$ \frac{\partial f}{\partial t} = \nabla \cdot \left[ \left( D_f \nabla f - \frac{1}{4} \mu_f f \left\{ (1 - \rho(u)) \left( \mathbf{S}_{ff} \nabla f + \mathbf{S}_{fm} \nabla m + \mathbf{S}_{fe} \nabla e \right. \right. \right. \\ \left. \left. + f \left( \frac{\partial}{\partial g} \mathbf{S}_{ff} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{ff} \nabla e \right) + m \left( \frac{\partial}{\partial g} \mathbf{S}_{fm} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{fm} \nabla e \right) \right. \right. \\ \left. + e \left( \frac{\partial}{\partial g} \mathbf{S}_{fe} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{fe} \nabla e \right) \right) - (\mathbf{S}_{ff} f + \mathbf{S}_{fm} m + \mathbf{S}_{fe} e)(w_g \nabla g \right. \\ \left. + w_f \nabla f + w_m \nabla m + w_e \nabla e) \right] - \lambda_f f + p_f(g) f(1 - \rho(u)), \qquad (2.22b) $$

$$ \frac{\partial m}{\partial t} = \nabla \cdot \left[ D_m \nabla m - \frac{1}{4} \mu_m m \left\{ (1 - \rho(u)) \left( \mathbf{S}_{mf} \nabla f + \mathbf{S}_{mm} \nabla m + \mathbf{S}_{me} \nabla e \right. \right. \right. \\ \left. \left. + f \left( \frac{\partial}{\partial g} \mathbf{S}_{mf} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{mf} \nabla e \right) + m \left( \frac{\partial}{\partial g} \mathbf{S}_{mm} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{mm} \nabla e \right) \right. \right. \\ \left. + e \left( \frac{\partial}{\partial g} \mathbf{S}_{me} \nabla g + \frac{\partial}{\partial e} \mathbf{S}_{me} \nabla e \right) \right) - (\mathbf{S}_{mm} m + \mathbf{S}_{mf} f + \mathbf{S}_{me} e)(w_g \nabla g \right. \\ \left. + w_f \nabla f + w_m \nabla m + w_e \nabla e) \right] - \lambda_m m + p_m(g) m(1 - \rho(u)), \qquad (2.22c) $$

$$ \frac{\partial e}{\partial t} = -e(\alpha_f f + \alpha_m m) + p_e(f)e(1 - \rho(u)). \qquad (2.22d) $$
### 3. Numerical approach and simulations

Throughout this section, we investigate numerically the macrophage-ECM hypothesis mentioned above (see cases I and II). We do this for both the original non-local system (2.1) and the local versions of it: systems (2.17) and (2.22).

#### 3.1. The finite element discretisation

Denoting by $\Omega$, the 2D spatial domain of interest and $I = [0, T]$ the temporal domain, we seek to approximate here the solution $g, f, m, e$ over $H^1(\Omega; I)$. For a concise presentation of the discretisation of Eq (2.1), let's denote by $u := (g, f, m, e)^T$ and the right-hand side operator by $G : (H^1(\Omega; I))^4 \to \mathbb{R}^4$. Thus, our dynamics re-casts as:

$$ egin{cases} 
\frac{\partial u}{\partial t} = G(u), \\
u(\mathbf{x}, 0) = u_0, \\
\left. \frac{\partial u}{\partial n} \right|_{\partial\Omega} = 0.
\end{cases} 
\quad (3.1) $$

For the spatial discretisation, we use the finite element method (FEM). Thus, as usual, we multiply the equation by a test function $v \in \mathcal{D}(\Omega)$ (the usual space of test functions, which coincides with $C_0^\infty(\Omega)$ as family of functions only) and integrate over $\Omega$. Therefore, in weak formulation, our problem re-casts as follows:

find $u \in (H^1(\Omega; I))^4$ such that:

$$ egin{cases} 
\displaystyle \int_{\Omega} \frac{\partial u}{\partial t} v \, d\mathbf{x} = \int_{\Omega} G(u) v \, d\mathbf{x}, & \forall v \in C_0^\infty(\Omega) \\
u(\mathbf{x}, 0) = u_0, \\
\left. \frac{\partial u}{\partial n} \right|_{\partial\Omega} = 0,
\end{cases} 
\quad (3.2) $$

where $n$ is the usual normal to $\Omega$. Thus, explicitly, the first equation in (3.2) becomes

$$ \int_{\Omega} \frac{\partial g}{\partial t} v_g \, d\mathbf{x} = - \int_{\Omega} D_g \nabla g \cdot \nabla v_g \, d\mathbf{x} - \int_{\Omega} \lambda_g g v_g \, d\mathbf{x} + \int_{\Omega} W(f, m) v_g \, d\mathbf{x} \\ + \int_{\partial\Omega} D_g \frac{\partial g}{\partial n} v_g \, dS, \quad (3.3a) $$

$$ \int_{\Omega} \frac{\partial f}{\partial t} v_f \, d\mathbf{x} = - \int_{\Omega} D_f \nabla f \cdot \nabla v_f \, d\mathbf{x} + \int_{\Omega} \nabla \cdot (f A_f [g, f, m, e]) v_f \, d\mathbf{x} - \int_{\Omega} \lambda_f f v_f \, d\mathbf{x} \\ + \int_{\Omega} p_f(g) f (1 - \rho(u)) v_f \, d\mathbf{x} + \int_{\partial\Omega} D_f \frac{\partial f}{\partial n} v_f \, dS, \quad (3.3b) $$

$$ \int_{\Omega} \frac{\partial m}{\partial t} v_m \, d\mathbf{x} = - \int_{\Omega} D_m \nabla m \cdot \nabla v_m \, d\mathbf{x} + \int_{\Omega} \nabla \cdot (m A_m [g, f, m, e]) v_m \, d\mathbf{x} $$
$$ -\int_{\Omega} \lambda_m m v_m \, dx + \int_{\Omega} p_m(g)m(1-\rho(u))v_m \, dx + \int_{\partial\Omega} D_m \frac{\partial m}{\partial n} v_m \, dS \quad (3.3c) $$

$$ \int_{\Omega} \frac{\partial e}{\partial t} v_e \, dx = - \int_{\Omega} e(\alpha_f f + \alpha_m m) v_e \, dx + \int_{\Omega} e p_e (1 - \rho(u)) v_e \, dx. \quad (3.3d) $$

Thus, since for each component $u_i$ of $u$, $i = 1, \dots, 4$, a representation in terms of P2-basis functions $\{\psi_{\tau}(\cdot)\}_{\tau=0,l}$ corresponding to the rectangular grid (consisting of $l+1$ uniformly distributed nodes) is given as $\tilde{u}(\mathbf{x}, t) = \sum_{\tau=0}^{l} c_{\tau}^{u_i}(t) \psi_i(\mathbf{x})$, we therefore obtain

$$ \frac{\partial u_i(\mathbf{x}, t)}{\partial t} = \sum_{\tau=0}^{l} \frac{\partial c_{\tau}^{u_i}(t)}{\partial t} \psi_{\tau}(\mathbf{x}). \quad (3.4) $$

For the time discretization, we use a standard backward Euler scheme, and thus from Eq (3.2) we obtain

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \frac{u^{N+1} - u^N}{\Delta t} v \, dx = \int_{\Omega} G(u^{N+1}) v \, dx, \quad \forall v \in \mathcal{D}(\Omega) \\[1.5ex] \displaystyle u(\mathbf{x}, 0) = u_0, \\[1.5ex] \displaystyle \frac{\partial u^{N+1}}{\partial n} = 0, \end{array} \right. \quad (3.5) $$

where $u^N$ is the approximation of $u(N\Delta t)$, with the uniform time discretisation step $\Delta t := \frac{T}{N_{max}}$, with $N = 0, \dots, N_{max}$ representing the time indices. Further, in order to finally write our discretised dynamics in the standard form of a linear system of equations, we first proceed with taking the $L^2$-scalar product with respect to each basis function $\psi_j$, $j = 0, \dots, l$, in both sides of equations in (3.3), and then we apply the time discretisation outlined in (3.5). This results in a linear system associated with the fully discretised model, as detailed in Appendix 4.

## 3.2. Numerical simulations

**Numerical implementation.** For the numerical simulations, we consider the finite element method implemented in FENiCS, an open-source computing platform. In particular, here we use a triangular mesh with P2 elements (i.e., on every single side (edge) of the triangles in the mesh the solution is approximated by a quadratic that is sampled at the following three points: the two vertices and the middle point); for precise details, see [52]. For the details of the weak form discretisation, see sub-section 3.1 and Appendix C. For the time-marching using the backward Euler, see Appendix D. We run the simulations on a time interval $[0, T]$ with $T = 100$ and a step size $\Delta t = 0.2$.

Finally, the parameter values used for the numerical simulations are listed in Table 1, together with their description. We assume homogeneous Neumann boundary conditions on all sides of the spatial domain $\Omega$.

**Initial conditions.** To properly model the healing of a wound (i.e., a decrease in the normal density of cells due to a cut in the tissue), we choose the following initial conditions on a square domain $\Omega = [-1, 1]^2$:

$$ g(\mathbf{x}, 0) = 0.1, \quad (3.6) $$

$$ f(\mathbf{x}, 0) = 0.4[(0.5 + 0.5 \tanh(20x_1 - 3)) + (0.5 + 0.5 \tanh(-20x_1 - 3))], \quad (3.7) $$
$$m(\mathbf{x}, 0) = 0.1[(0.5 + 0.5 \tanh(20x_1 - 3)) + (0.5 + 0.5 \tanh(-20x_1 - 3))], \quad (3.8)$$

$$e(\mathbf{x}, 0) = 1.0[(0.5 + 0.5 \tanh(20x_1 - 3)) + (0.5 + 0.5 \tanh(-20x_1 - 3))], \quad (3.9)$$

where $\mathbf{x} = (x_1, x_2)$. Thus the wound is assumed to be parallel to the $x_2$-axis; see also Figure 2.

<table><caption>**Table 1.** Summary of dimensionless model parameters, together with the baseline values used for the numerical simulations corresponding to limit local models.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Description</th><th>Reference</th></tr></thead><tbody><tr><td>D<sub>g</sub></td><td>0.0035</td><td>Diffusion coeff. for growth-factor population</td><td>[53]</td></tr><tr><td>D<sub>f</sub></td><td>0.0008</td><td>Diffusion coeff. for fibroblast population</td><td>[53]</td></tr><tr><td>D<sub>m</sub></td><td>0.0008</td><td>Diffusion coeff. for macrophage population</td><td>[53]</td></tr><tr><td>λ<sub>g</sub></td><td>0.2</td><td>Decay rate of growth-factor population</td><td>Estimated</td></tr><tr><td>λ<sub>f</sub></td><td>0.025</td><td>Apoptotic rate of fibroblast population</td><td>Estimated</td></tr><tr><td>λ<sub>m</sub></td><td>0.025</td><td>Apoptotic rate of macrophages population</td><td>Estimated</td></tr><tr><td>p<sup>f</sup><sub>g</sub></td><td>0.2</td><td>Secretion rate of growth-factor by fibroblasts</td><td>Estimated</td></tr><tr><td>p<sup>m</sup><sub>g</sub></td><td>0.2</td><td>Secretion rate of growth-factor by macrophages</td><td>Estimated</td></tr><tr><td>p<sub>f</sub>(g)</td><td>5.0g</td><td>Proliferation rate of fibroblasts population depending on the growth factor</td><td>Estimated</td></tr><tr><td>p<sub>m</sub>(g)</td><td>5.0g</td><td>Proliferation rate of macrophages population depending on the density of the growth factor</td><td>Estimated</td></tr><tr><td>α<sub>f</sub></td><td>0.015</td><td>Degradation rate of ECM by fibroblasts</td><td>[53]</td></tr><tr><td>α<sub>m</sub></td><td>0.015</td><td>Degradation rate of ECM by macrophages</td><td>[53]</td></tr><tr><td>p<sub>e</sub>(f)</td><td>5.0 f</td><td>Remodelling rate of ECM population</td><td>Estimated</td></tr><tr><td>w<sub>g</sub></td><td>1</td><td>Fraction of physical space occupied by growth factor</td><td>[54]</td></tr><tr><td>w<sub>f</sub></td><td>1</td><td>Fraction of physical space occupied by fibroblasts</td><td>[54]</td></tr><tr><td>w<sub>m</sub></td><td>1</td><td>Fraction of physical space occupied by macrophages</td><td>[54]</td></tr><tr><td>w<sub>e</sub></td><td>1</td><td>Fraction of physical space occupied by ECM</td><td>[54]</td></tr><tr><td>S<sup>max</sup><sub>ff</sub></td><td>0.2</td><td>Maximum strength of fibroblast-fibroblast adhesive junction</td><td>Estimated</td></tr><tr><td>S<sup>max</sup><sub>fm</sub></td><td>0.1</td><td>Maximum strength of fibroblast-macrophages adhesive junction</td><td>[51]</td></tr><tr><td>S<sup>max</sup><sub>mf</sub></td><td>0.1</td><td>Maximum strength of macrophages-fibroblast adhesive junction</td><td>[51]</td></tr><tr><td>S<sup>max</sup><sub>mm</sub></td><td>0.2</td><td>Maximum strength of macrophages-macrophages adhesive junction</td><td>Estimated</td></tr><tr><td>S<sup>max</sup><sub>fe</sub></td><td>0.1</td><td>Maximum strength of fibroblast-ECM adhesive junction</td><td>[51]</td></tr><tr><td>S<sup>max</sup><sub>me</sub></td><td>1.0</td><td>Maximum strength of macrophages-ECM adhesive junction</td><td>[51]</td></tr><tr><td>μ<sub>f</sub></td><td>0.08</td><td>Haptotactic rate of the fibroblasts</td><td>Estimated</td></tr><tr><td>μ<sub>m</sub></td><td>0.08</td><td>Haptotactic rate of the macrophages</td><td>Estimated</td></tr><tr><td>R</td><td>0.1</td><td>Sensing radius for the non-local interaction</td><td>[51]</td></tr></tbody></table>
Note that in Appendix F we show also numerical simulations for the baseline parameters when we choose different initial conditions: a circular initial wound (see Eq (4.13d) and first column in Figure A3) and a linear but irregular initial wound (see Eq (4.12d) and first column in Figure A4).

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Initial conditions (see Eqs (3.6)–(3.9)) for: (a) growth factor, *g*; (b) fibroblasts, *f*; (c) macrophages, *m* and (d) ECM, *e*. We consider a square domain Ω = [-1, 1]².</figcaption></figure>

3.2.1. Numerical simulations for the local models: Gaussian kernel

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Numerical simulations of the *local model* (2.17), obtained as a limit from the non-local model with a Gaussian kernel. The rows correspond to the spatial distribution of growth factor (*g*), fibroblast (*f*), macrophages (*m*), and ECM (*e*) at times *t* = 2, *t* = 20, *t* = 100 and *t* = 200 (shown on the columns). The parameter values are listed in Table 1.</figcaption></figure>
We start the numerical simulations by focusing first on local models, and in particular on Eq (2.17). Here, the motility of cells is only due to diffusion as the cellular flux vanishes locally (see the detailed calculations in Appendix 4). Figure 3 shows the spatio-temporal evolution of growth factor *g* (first row), fibroblasts *f* (second row), macrophages *m* (third row), and ECM *e* (fourth row) at three different time points: *t* = 2, *t* = 20, *t* = 100 and *t* = 200.

We see in Figure 3 that due to diffusion, there is an invasion of fibroblasts and macrophages into the wound between their initial states (*t* = 0, see Figure 2) and time *t* = 2, which helps also remodel the ECM. Note that at time *t* = 20, the fibroblast and macrophages have completely invaded the wound and now started decaying. For the parameter values used here (and listed in Table 1), the ECM returns to its maximum level, while the fibroblast and macrophage tissue levels decay compared to their initial baseline level (at *t* = 0).

3.2.2. Numerical simulations for the local models: cone-shaped kernel

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Numerical simulations of the *local model* (2.22) with $S_{me} = 0$ obtained as a limit from the non-local model Case I with a cone-shaped kernel. The rows correspond to the spatial distribution of growth factor (*g*), fibroblasts (*f*), macrophages (*m*), and ECM (*e*) at time points *t* = 2, *t* = 20, and *t* = 100 (shown on the columns). The parameter values are listed in Table 1.</figcaption></figure>
**No macrophage-ECM adhesion (Case I: S<sub>me</sub> = 0).** First, we investigate numerically the behaviour of the local model (2.22) under the assumption that the macrophages cannot adhere to the extracellular matrix (ECM), i.e., case I above with S<sub>me</sub> = 0. As before, in Figure 4, we plot the spatio-temporal evolution of growth factor *g* (first row), fibroblasts *f* (second row), macrophages *m* (third row), and ECM *e* (fourth row). In contrast to Figure 3 where the fibroblasts were smoothly invading the wound area, in Figure 4, we see that at *t* = 2, the fibroblasts invade the wound area in oscillatory waves. Then, at *t* = 20, there is a very large density of fibroblasts built-up in the middle of the wound (*max(f)* ≈ 0.5) and a slightly lower fibroblast density at the wound edges. This very high peak of fibroblasts then starts decreasing and at *t* = 100, the wound is almost healed and the fibroblasts are almost homogeneously spread over the whole tissue at very low densities. For the parameter values used here (and listed in Table 1), the ECM returns to its maximum level, while the fibroblast tissue level decays compared to its initial baseline level (at *t* = 0).

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Numerical simulations of the *local model* (2.22), obtained as a limit from the non-local model in Case II ($S_{me} > 0$), with a cone-shaped kernel. The rows correspond to the spatial distribution of growth factor (*g*), fibroblast (*f*), macrophages (*m*), and ECM (*e*) at time points *t* = 2, *t* = 20, and *t* = 100. The parameter values are listed in Table 1.</figcaption></figure>

**Including macrophage-ECM adhesion (Case II: S<sub>me</sub> > 0).** Next, we investigate numerically the case where macrophages could adhere to the ECM. In contrast to the previous figure (i.e., Figure 4), we see in Figure 5 that at *t* = 2, there is an invasion of macrophages into the wound (in an oscillatory manner) as a result of cell-ECM adhesion, which helps remodel also the ECM. Another interesting macrophage behaviour can be observed at time *t* = 20, when there are fewer macrophages inside the
wound (slightly darker colour) compared to the edge of the wound and the rest of the tissue (see row 3 column 2 of Figure 5). This is in contrast to Figures 3 and 4. For the parameter values used here (and listed in Table 1), the ECM in the wound region remodels faster and it reaches its maximum level at $t = 100$. This is in contrast to Figures 4 and 3, where ECM remodelling is slightly slower.

Finally, in Figure 6, we compare the solutions of these two hypotheses (i.e., Figure 4 with $S_{me} = 0$, and Figure 5 with $S_{me} > 0$) at the spatial point $x_2 = 0$ and times $t = 2, 20, 40, 100$. It is clear that the ECM remodels faster when macrophages are allowed to adhere to the ECM.

**Time-evolution of fibroblasts & macrophages: comparison between different local models.** To see more clearly the differences between various local models in terms of fibroblasts and macrophages, in Figure 7 we plot the densities of these two cell populations in the center of the wound (i.e., at point $x_1 = 0$; blue curve) and just outside the wound (i.e., at point $x_1 = 0.5$; red curve) for time $t \in \{0, 50\}$. Here we show cell dynamics for: (a),(b) local model (2.17); (c),(d) local model (2.22) for case I ($S_{me} = 0$); (e),(f) local model (2.22) for case II ($S_{me} > 0$). In Figure 7(a),(b) we observe a similar dynamics for fibroblasts and macrophages: for $t > 5$ there are more cells in the wound region compared to surrounding tissue. This is in contrast to Figure 7(c),(d) where the fibroblasts population invades very quickly the wound ($t > 3$) and reaches high levels very quickly, while the macrophage population in the wound region does not exceed that of the surrounding tissue until about $t \in \{5, 10\}$. Finally, in Figure 7(e),(f), the fibroblast population in the wound region exceeds that of the surrounding tissues as early as $t > 4$, and it reaches its maximum peak around $t = 25$. The macrophage population in the wound region exceeds the level of macrophages in the surrounding tissue at $t = 3$, but then quickly decreases below the level of macrophages in the surrounding tissue around $t = 15$.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** 2-D plots comparing the solutions of the *local models* obtained via the localization of cone-shaped kernel for Cases I and II. We show the spatial distribution of variables at: (a) $t = 2$; (b) $t = 20$; (c) $t = 40$ and (d) $t = 100$. The parameter values are listed in Table 1.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The evolution of fibroblasts (left column) and macrophages (right column) inside the wound (at $x_1 = 0$) and in the neighbouring tissue (at $x_1 = 0.5$) for the *local models*. Here we compare the effect of various kernels on the localised models: (a),(b) Gaussian kernel; (c),(d) cone-shaped kernel (Case I: $S_{me} = 0$); and (e),(f) cone-shaped kernel (Case II: $S_{me} > 0$), from $t = 0$ to $t = 50$. The parameter values are listed in Table 1.</figcaption></figure>

### 3.2.3. Numerical simulations for the non-local models

Next, we investigate numerically the non-local model (2.1), while focusing again mostly on Case II (where macrophages-ECM adhesion is present). For the parameter values used here (and listed in Table 1), the ECM returns to its maximum level, while the fibroblasts level decays compared to its initial baseline level (at $t = 0$).

Figure 8 shows the spatio-temporal evolution of growth factor *g* (first row), fibroblasts *f* (second row), macrophages *m* (third row), and ECM *e* (fourth row) for the non-local model with cone-shaped kernel, when macrophage-ECM adhesion is allowed (Case II). We see here that as early as $t > 2$, fibroblasts and macrophages invade the wound as a result of cell diffusion and adhesion, which helps also remodel the ECM. At time $t = 20$ there are many more fibroblasts inside the wound compared to wound margins where there are low numbers of fibroblasts.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Numerical simulations of the *non-local model* (2.1) with the cone-shaped kernel (2.6), for case II (macrophage-ECM adhesion present). The rows correspond to the spatial distribution of growth factors (*g*), fibroblasts (*f*), macrophages (*m*), and ECM (*e*) at times *t* = 2, *t* = 20, and *t* = 100 (on the columns). The parameter values are listed in Table 1.</figcaption></figure>

In Figure 9 we compare cases I ($S_{me} = 0$) and II ($S_{me} > 0$) for the non-local model (2.1) with (a)–(d) a cone-shaped kernel, and (a')–(d') a Gaussian kernel. We see that there are no significant differences between cases I. and II. Only at *t* = 100 we can observe a slightly faster ECM remodelling for case I compared to case II (although this difference is barely noticeable, probably because $S_{me}$ is not large enough). Moreover, there are no differences between the results obtained with cone-shaped kernel (top panels) and those obtained with the Gaussian kernel (bottom panels).

To confirm this similarity in the results, in Figure 10 we show a log-log plot of the $L_2$ norm of the difference between the solution obtained with the Gaussian kernel and the solution obtained the cone-shaped kernel (for Case II: $S_{me} > 0$) at each time step, i.e., from *t* = 0 to *t* = 100. In this figure, we see that $L_2$ norm of the difference in the solutions is mostly between $10^{-7} - 10^{-6}$, confirming that the two solutions are equal when approximated to 6 decimal places.

Because of these similarities in the solutions obtained with Gaussian and cone-shaped kernels, throughout the rest of this study we consider only cone-shaped kernels.

### 3.2.4. Local vs non-local models

Next, we highlight some differences between the local model (2.22) and non-local model (2.1), when we fix all parameters. To this end, we focus only on Case II (macrophage-ECM adhesion present) and on cone-shaped kernels.

In Figure 11(a) we see that as early as *t* = 2 fibroblasts and macrophages invade the wound region via random movement (diffusion) and directed movement (haptotaxis), and this invasion is charac-
Cone-shaped Kernel:

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

Gaussian Kernel:

<figure><img src="image_12.png" /><figcaption>**Figure 9.** 2-D plots comparing the solution of the *non-local model* (2.1) with the cone-shaped kernel (sub-panels (a)–(d)), and the Gaussian kernel (sub-panels (a')–(d')) for Case I ($S_{me} = 0$) and Case II ($S_{me} > 0$). We show the spatial distribution of variables at: (a) $t = 2$; (b) $t = 20$; (c) $t = 40$ and (d) $t = 100$. The parameter values are listed in Table 1.</figcaption></figure>

terised by oscillations in ECM and fibroblasts densities at the wound margin. These oscillations seem
to have slightly higher amplitudes for the local model compared to the non-local one. As time in-
creases, these oscillations die out and the ECM repairs. However, there is an interesting difference
between the local and non-local models:

* At $t = 40$: for the local model (2.22) the ECM is still very low at the initial wound point ($x_1 = 0$), while for the non-local model (2.1) the ECM has much higher densities at $x_1 = 0$.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Log-log representation of the $L_2$ norm of the difference between the solution of the *non-local model* (2.1) with Gaussian kernel, and the solution of the same non-local model with cone-shaped kernel. Here we focus only on Case II (macrophage-ECM adhesion present). The $L_2$ norm is calculated at every time point: $t = 0$ to $t = 100$.</figcaption></figure>

* At $t = 100$: for the local model the ECM has almost remodelled completely (especially at $x_1 = 0$), while for the non-local model the ECM is still remodelling and filling up the wound.

We suspect that this faster ECM remodelling (and wound healing) observed with the local model is the result of higher oscillation amplitudes in the ECM densities at the wound margins (particularly evident at $t = 40$).

<figure><img src="image_14.png" /><figcaption>**Figure 11.** 2-D plots comparing the solution of the *local model* to the *non-local model* for case II (i.e., where macrophage-ECM adhesion is present) using the cone-shaped kernel at the time (a) $t = 2$; (b) $t = 20$; (c) $t = 40$ and (d) $t = 100$. The parameter values are listed in Table 1.</figcaption></figure>

Note that in addition to the above numerical investigation of the impact of cell sensing radius ($R = 0$ for local case, $R = 0.1$ for non-local case) on the healing of the wound, in Appendix G we have included also a set of simulations showing model dynamics for $R = 0.08$, $R = 0.1$ (baseline, investigated also above), and $R = 0.13$.
### 3.3. Assumption of no cell death due to overcrowding

All simulations in the previous section showed normal wound healing. To highlight one possible biological mechanism that could lead to abnormal wound healing characterised by raised scars (as seen in hypertrophic or keloid scars), in the following we investigate numerically the hypothetical case where cells do not die due to overcrowding. In fact, in [55] it was observed experimentally that the fibroblasts in the center of the keloid lesion had reduced doubling time and also lower death rates, and thus reached higher cell densities compared to the saturated-like densities of normal fibroblasts. Because of these observed higher fibroblasts densities with reduced death rates, in the following we replace the classical logistic proliferation terms appearing in Eqs (2.1b) and (2.1c) with the following truncated-logistic terms that ignore cell death at higher densities (see also Appendix 4):

$$p_f f(1 - \rho(u))^+ \quad \text{and} \quad p_m m(1 - \rho(u))^+. \qquad (3.10)$$

Figures 12–14 show model dynamics for this particular case in the context of local and non-local interactions.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Numerical simulations of the *local model* (2.22) case II, with a cone-shaped kernel, for the hypothetical case where cells do not die because of overcrowding (see Eq (3.10)). The rows correspond to the spatial distribution of the growth factor (*g*), fibroblasts (*f*), macrophages (*m*), and ECM (*e*), at times *t* = 2, *t* = 20, and *t* = 100. The parameter values are listed in Table 1.</figcaption></figure>

* In Figure 12, we show the dynamics of the *local model* (2.22) (for case II: $S_{me} > 0$), with these new proliferation rates when we localize the cone-shaped kernel. We note that the dynamics is similar to the dynamics in Figure 5: a normal wound healing.
* In Figure 13, we show the dynamics of the non-local model (2.1) (for case II: $S_{me} > 0$) with cone-shaped kernel. In contrast to Figure 12, here we consider also $e_c > 0$ describing the possibility of cells moving down ECM gradients in the first stages of the wound healing. We see that the concentration of the growth factor and the density of fibroblasts grow significantly at $t = 100$, eventually leading to the blow-up of the numerical code. The growth in the fibroblasts population is not matched by a similar growth in macrophages; this could be an indirect result of the non-linear interactions and the asymmetry in the fibroblasts-ECM and macrophage-ECM interactions (where only fibroblasts are assumed to contribute to ECM remodelling).

Since the non-local model (2.1) with cone-shaped kernel, truncated logistic cell growth and $e_c > 0$ in the adhesion function leads to very high fibroblasts densities (as seen in Figure 13), next we investigate whether this fibroblasts dynamics holds also for a Gaussian kernel. In Figure 14 we see that while the non-local model with Gaussian kernel exhibits the same overgrowth of fibroblasts (for the same parameter values as in Figure 13), the corresponding local model does not show fibroblasts overgrowth.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Numerical simulations of the *non-local model* (2.1) with *truncated logistic term for both cells and the ECM*. We consider case II with a cone-shaped kernel. The rows correspond to the spatial distribution of growth factor (*g*), fibroblast (*f*), macrophages (*m*), and ECM (*e*) at time points *t* = 2, *t* = 20, and *t* = 100. The parameter values are listed in Table 1 with the following adjustments: $\mu_f = \mu_m = 10.0$, $p_f(g) = 20g$, $\lambda_f = 0.0000025$, and $e_c = 0.9$.</figcaption></figure>

### 3.4. Primary vs. secondary wound healing

As mentioned in the Introduction section 1, wound healing can occur by primary intention (when the woundheals as the wound margins are coming together; as it is the case of surgical incisions, skin grafts, or flap closures) or by secondary intention (when the wound is very large and itheals from the bottom up as the granulation tissue is formed and fills in the wound) [56]. In Figure 8, we observed wound healing by secondary intention, as the ECM was remodelled from the bottom up. In contrast,
<figure><img src="image_17.png" /><figcaption>**Figure 14.** 2-D plots comparing the solution of the *local model* to the *non-local model* for case II (i.e., macrophage-ECM adhesion is present: $S_{me} > 0$) using a Gaussian kernel and a *truncated logistic term for both cells and the ECM* (a) $t = 2$; (b) $t = 20$; (c) $t = 40$ and (d) $t = 100$. The parameter values are listed in Table 1 with the following adjustments: $\mu_f = \mu_m = 10.0$, $p_f(g) = 20g$, $\lambda_f = 0.0000025$, and $e_c = 0.9$</figcaption></figure>

in Figure 15 we show wound healing by primary intention as the wound closes from the sides; see the ECM progression from $t = 2$ to $t = 20$ and finally $t = 100$. To obtain this dynamics, we reduced the diffusion rate of fibroblasts.

**4. Summary and discussion**

In this study, we developed a new mathematical model to describe some simple interactions between fibroblasts, macrophages, ECM and a growth factor in the context of wound healing. Due to the non-local aspects of the cell-cell and cell-ECM interactions (which are the result of various factors, from adhesive forces [45, 46, 51, 53] to non-conventional cell protrusions that allow long-distance cell-cell interactions [23, 24]), we started with a non-local model that considers a non-local flux generated by these bio-mechanical attractive/adhesive/repulsive interactions. However, since the spatial range over which cells perceive other neighbouring cells is not always very clear, we also considered localised versions of the original non-local model by assuming that cell perception radius ($R$) approaches zero. We showed that the type of kernel that describes the non-local interactions has a significant impact on the resulting local models.

More precisely, Gaussian kernels led to local reaction-diffusion models, while cone-shaped kernels led to local reaction-advection-diffusion models. In addition, since the published literature is not very clear about the adhesion of macrophages to ECM, throughout this study we investigated numerically two possible cases: (I.) no macrophage-ECM adhesion ($S_{me} = 0$); and (II.) macrophage-ECM adhesion present ($S_{me} > 0$). Finally, since the wound starts to heal with the formation of the fibrin mesh that acts as an early ECM, to describe the initial movement of cells into the wound we chose two different adhesion functions.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Numerical simulations of the *non-local model* (2.1) with the *classical logistic term for both cells and the ECM*. We consider case II with a cone-shaped kernel. The rows correspond to the spatial distribution of growth factor (g), fibroblast (f), macrophages (m), and ECM at times points *t* = 2, *t* = 20, and *t* = 100, corresponding to the columns respectively. The parameter values are listed in Table 1 with *D*<sub>*f*</sub> = 0.000008.</figcaption></figure>

Numerical simulations were performed for all models (local and non-local, with/without macrophage-ECM adhesive interactions), to illustrate some of the behaviours exhibited by these systems. The results showed that ECM remodelling depends on the type of kernels considered and on the local vs. non-local adhesive interactions. More precisely, the results showed that ECM remodelling is slower in the non-local models compared to the local models, ECM remodelling is slower in the presence of macrophage-ECM adhesion compared to when macrophage-ECM is not present in both the local and non-local models. The results also showed that a reduction in the diffusion of fibroblasts in the non-local model may lead to wound healing by primary intention (see Figure 15), while a higher fibroblast diffusion is associated with healing by secondary intention (see Figure 8).

Another interesting result is the abnormal healing observed in the non-local model (Figures 13 and 14) due to the uncontrolled proliferation of fibroblasts and production of growth factor, in a specific situation: when (i) cells do not die due to overcrowding, and (ii) cells move down-gradient at the beginning of the healing process, onto the fibrin mesh that fills in the wound (process modelled by the introduction of an ECM threshold $e_c > 0$). Since the fibroblasts density is above the threshold $f = 1$ throughout a spatial region that exceeds the original wound region, we can say that this case corresponds to keloid scar formation.

Also interesting is the transient oscillatory invasive patterns of fibroblasts into the wound gap (see row 2 column 1 of Figures 4 and 5). We note that oscillations of cell density fluctuations have been observed in experiments on epithelial tissues [57], where they have been suggested to be the result of
cell-cell adhesion.

To conclude, we emphasise that all these results are purely theoretical, showing the possible model behaviours as we vary randomly the non-dimensional parameters. Comparison with experimental data is necessary to *quantitatively* investigate biologically-relevant normal and abnormal wound healing patterns. Nevertheless, the simulation results in this study can describe *qualitatively* many of the observed normal and abnormal wound healing processes: from the spatial collective oscillators in the cell density (Figures 4 and 5) observed experimentally in the proliferation stage as the wound tissue heals [57] (note also that some clinical studies refer to the proliferation stage as the granulation tissue formation stage due to the granular appearance of the tissue [58], that can also lead to an oscillatory wave-like appearance of the tissue); to the higher ECM densities sometimes observed temporary at the site of the wound (see Figures A3 and A4 in Appendix F) which then subside in the absence of immune cells [59] as the tissue heals normally; the higher fibroblast densities (see Figure 13) observed in abnormal wound healing that gives rise to hypertrophic or keloid scars; the healing by primary intention (Figure 15) or by secondary intention (Figure 8). While the general shapes of these numerically-simulated wounds can be qualitatively compared with the shapes of actual clinical wounds, we cannot claim that our simple mathematical model (that does not consider, for example, the heterogeneity of macrophages or fibroblasts inside the normal/abnormal wounds [60, 61]) captures accurately all aspects involved in a clinical wound healing.

**Future work.** First, as mentioned above, we need to estimate model parameters using real 2D experimental data (using inverse problem approaches, as in [62]). In addition, we also need to estimate from the data the right kernels (for the spatial ranges of the non-local interactions), the cell-cell and cell-ECM adhesion functions, and proliferation laws. However, at this moment, we do not have such detailed immunological and bio-mechanical data (and we could not find it in the published literature).

Second, this theoretical study generated a few more theoretical questions that will be investigated in the future. More precisely, the numerical simulations in Figure 3 (last two columns) suggested that the solution likely approaches some spatially homogeneous steady states with various magnitudes. Hence, it is normal to investigate the existence and linear stability of such solutions characterised by the spatially-uniform spread of cell/molecule densities across all domain, to gain a better understanding of the parameters (and biological mechanisms) that are behind the various behaviours. Such an analysis would also allow us to possibly identify the biological mechanisms that could contribute to abnormal wound healing (via instability of the analytically-identified spatially-homogeneous states and stability of the spatially-heterogeneous states). Moreover, some numerical simulations (not shown here) showed numerical instabilities in the solutions, followed by numerical blow-up. This is an expected outcome for advection-dominated models discretised using FEM. In the future is important to stabilise these numerical schemes to ensure more accurate results (especially if we want to compare with experimental data).

Third but not last, we will investigate the well-posedness of local and non-local models introduced in this study, which will also ensure the local regularity of the solutions. This analysis will complement the current numerical investigation, to elucidate the mechanisms behind the numerically-simulated abnormal wound healing behaviour (characterised by excessive growth in the densities of some model components; see Figure 13).
**Use of AI tools declaration**

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

**Acknowledgments**

RE and OA acknowledge funding from the French Agence Nationale de la Recherche (ANR) grant number ANR-21-CE45-0025-01; SPAB and SU acknowledge funding from the Luxembourg National Research Fund (FNR) grant number INTER/ANR/21/16399490; G.R. acknowledges funding from the ANR grant number ANR-21-CE45-0025-03.

**Conflict of interest**

The authors declare there is no conflict of interest.

**References**

1. R. F. Diegelmann, R. F. Evans, Wound healing: an overview of acute, fibrotic and delayed healing, *Front. Biosci.*, **9** (2004), 283–289. https://doi.org/10.2741/1184

2. D. Zuo, D. He, H. Yang, K. H. Stephane Avril, A thermodynamic framework for unified continuum models for the healing of damaged soft biological tissue, *J. Mech. Phys. Solids*, **158** (2022), 104662. https://doi.org/10.1016/j.jmps.2021.104662

3. S. Enoch, D. J. Leaper, Basic science of wound healing, *Surgery*, **23** (2005), 37–42. https://doi.org/10.1383/surg.23.2.37.60352

4. G. Gurtner, V. W. Wong, Wound Healing: Normal and Abnormal, Philadelphia, PA, (2014), 13–19.

5. G. C. Limandjaja, L. J. van den Broek, T. Waaijman, M. Breetveld, S. Monstrey, R. J. Scheper, et al., Reconstructed human keloid models show heterogeneity within keloid scars, *Arch. Dermatol. Res.*, **310** (2018), 815–826. https://doi.org/10.1007/s00403-018-1873-1

6. G. C. Limandjaja, F. B. Niessen, R. J. Scheper, S. Gibbs, Hypertrophic scars and keloids: Overview of the evidence and practical guide for differentiating between these abnormal scars, *Exp. Dermatol.*, **30** (2021), 146–161. https://doi.org/10.1111/exd.14121

7. S. GibbsPakyari, S. GibbsFarrokhi, M. K. Maharlooei, A. Ghahary, Critical role of transforming growth factor beta in different phases of wound healing, *Adv. Wound Care*, **2** (2013), 215–224. https://doi.org/10.1089/wound.2012.0406

8. S. Sanjabi, L. A. Zenewicz, M. Kamanaka, R. A. Flavell, Anti-inflammatory and pro-inflammatory roles of TGF-β, IL-10, and IL-22 in immunity and autoimmunity, *Curr. Opin. Pharmacol.*, **250** (2009), 447–53. https://doi.org/10.1016/j.coph.2009.04.008

9. E. Comellas, T. C. Gasser, T. C. Bellomo, S. Oller, A homeostatic-driven turnover remodelling constitutive model for healing in soft tissues, *J. R. Soc. Interface*, **13** (2016), 20151081. https://doi.org/10.1098/rsif.2015.1081
10. M. Pakyari, A. Farrokhi, M. K. Maharlooei, A. Ghahary, Critical role of transforming growth factor beta in different phases of wound healing, *Adv. Wound Care, **2*** (2013), 215-224.
https://doi.org/10.1089/wound.2012.0406

11. J. A. Flegg, J. A. Menon, P. K. Maini, D. L. S. McElwain, On the mathematical modeling of wound healing angiogenesis in skin as a reaction-transport process, *Front. Physiol., **6*** (2015), 1-17.
https://doi.org/10.3389/fphys.2015.00262

12. R. Eftimie, G. Rolin, O. Adebayo, S. Urcun, F. Chouly, S. P. A. Bordas, Modelling keloid dynamics: a brief review and new mathematical perspectives, Submitted, 2023.

13. J. A. Sherratt, J. D. Murray, Models of epidermal wound healing, *Proc. R. Soc. London, Ser. B, **241*** (1990), 29-36.
https://doi.org/10.1098/rspb.1990.0061

14. G. J. Pettet, H. M. Byrne, D. L. S. McElwain, J. Norbury, A model of wound-healing angiogenesis in soft tissue, *Math. Biosci., **136*** (1996), 35-63.
https://doi.org/10.1016/0025-5564(96)00044-2

15. G. Pettet, M. A. J. Chaplain, D. L. S. McElwain, H. M. Byrne, On the role of angiogenesis in wound healing, *Proc. R. Soc. London, Ser. B, **263*** (1996), 1487-1493.
https://doi.org/10.1098/rspb.1996.0217

16. E. A. Gaffney, K. Pugh, P. K. Maini, F. Arnold, Investigating a simple model of cutaneous wound healing angiogenesis, *J. Math. Biol., **45*** (2002), 337-374.
https://doi.org/10.1007/s002850200161

17. R. C. Schugart, A. Friedman, R. Zhao, C. K. Sen, Wound angiogenesis as a function of tissue oxygen tension: A mathematical model, *Proc. Natl. Acad. Sci., **105*** (2008), 26-28.
https://doi.org/10.1073/pnas.0711642105

18. M. Byrne, M. A. J. Chaplain, D. L. Evans, I. Hopkinson, Mathematical modelling of angiogenesis in wound healing: comparison of theory and experiment, *J. Theor. Med., **2*** (2000), 175-197.
https://doi.org/10.1080/10273660008833045

19. J. A. Flegg, H. M. Byrne, D. L. S. McElwain, Mathematical model of hyperbaric oxygen therapy applied to chronic diabetic wounds, *Bull. Math. Biol., **72*** (2010), 1867-1891.
https://doi.org/10.1007/s11538-010-9514-7

20. B. D. Cumming, D. L. S. McElwain, Z. Upton, A mathematical model of wound healing and subsequent scarring, *J. R. Soc. Interface, **7*** (2010), 19-34.
https://doi.org/10.1098/rsif.2008.0536

21. D. Zuo, S. Avril, H. Yang, S. J. Mousavi, K. Hackl, Y. He, 3D numerical simulation of soft tissue wound healing using constrained-mixture anisotropic hyperelasticity and gradient-enhanced damage mechanics, *J. R. Soc. Interface, **17*** (2020), 20190708.
https://doi.org/10.1098/rsif.2019.0708

22. Y. Kim, M. A. Stolarska, H. G. Othmer, A hybrid model for tumor spheroid growth in vitro I: Theoretical development and early results, *Math. Models Methods Appl. Sci., **17*** (2007), 1773-1798.
https://doi.org/10.1142/S0218202507002479

23. S. Caviglia, E. A. Ober, Non-conventional protrusions: the diversity of cell interactions at short and long distance, *Curr. Opin. Cell Biol., **54*** (2018), 106-113.
https://doi.org/10.1016/j.ceb.2018.05.013

24. D. S. Eom, Airinemes: thin cellular protrusions mediate long-distance signalling guided by macrophages, *Open Biol., **10*** (2020), 200039.
https://doi.org/10.1098/rsob.200039
25. C. Metzner, F. Hörsch, C. Mark, T. Czerwinski, A. Winterl, C. Voskens, Detecting long-lange interactions between migrating cells, *Sci. Rep.*, **11** (2021), 15031. https://doi.org/10.1038/s41598-021-94458-0

26. S. I. Despa, F. Despa, Diffusion model for growth factors—cell receptors interaction, *Biosystems*, **44** (1997), 59–68. https://doi.org/10.1016/S0303-2647(97)00047-6

27. F. Sefat, M. C. Denyer, M. Youseffi, Effects of different transforming growth factor beta-β isomers on wound closure of bone cell monolayers, *Cytokine*, **69** (2014), 75–86. https://doi.org/10.1016/j.cyto.2014.05.010

28. J. P. Andrews, J. Marttala, E. Macarak, J. Rosenbloom, J. Uitto, Keloids: The paradigm of skin fibrosis—pathomechanisms and treatment, *Matrix Biol.*, **51** (2016), 37–46. https://doi.org/10.1016/j.matbio.2016.01.013

29. P. Dicker, P. Pohjanpelto, P. Pettican, E. Rozengurt, Similarities between fibroblast-derived growth factor and platelet-derived growth factor, *Exp. Cell Res.*, **135** (1981), 221–227. https://doi.org/10.1016/0014-4827(81)90314-1

30. P. Sroobant, M. D. Waterfield, E. Rozengurt, Purification of fibroblast-derived growth factor, *Methods Enzymol.*, **147** (1987), 40–47. https://doi.org/10.1016/0076-6879(87)47097-3

31. A. Viola, F. Munari, R. Sánchez-Rodríguez, T. Scolaro, A. Castegna, The metabolic signature of macrophage responses, *Front. Immunol.*, **10** (2019), 1–16. https://doi.org/10.3389/fimmu.2019.00001

32. S. Huda, B. Weigelin, K. Wolf, K. V. Tretiakov, K. Polev, G. Wilk, et al., Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated *in vivo*, *Nat. Commun.*, **9** (2018), 4539. https://doi.org/10.1038/s41467-018-06563-w

33. R. J. Petrie, A. D. Doyle, K. M. Yamada, Random versus directionally persistent cell migration, *Nat. Rev. Mol. Cell Biol.*, **10** (2009), 538–549. https://doi.org/10.1038/nrm2729

34. M. V. Plikus, X. Wang, S. Sinha, E. Forte, S. M. Thompson, E. L. Herzog, et al., Fibroblasts: origins, definitions, and functions in health and disease, *Cell*, **184** (2021), 3852–3872. https://doi.org/10.1016/j.cell.2021.06.024

35. W. Jin, E. T. Shah, C. J. Penington, S. W. McCue, P. K. Maini, M. J. Simpson, et al., Logistic proliferation of cells in scratch assays is delayed, *Bull. Math. Biol.*, **79** (2017), 1028–1050. https://doi.org/10.1007/s11538-017-0267-4

36. S. Suveges, R. Eftimie, D. Trucu, Directionality of macrophages movement in tumour invasion: A multiscale moving-boundary approach, *Bull. Math. Biol.*, **82** (2020), 148. https://doi.org/10.1007/s11538-020-00819-7

37. G. C. Limandjaja, F. B. Niessen, R. J. Scheper, S. Gibbs, The keloid disorder: Heterogeneity, histopathology, mechanisms and models, *Front. Cell Dev. Biol.*, **8** (2020), 360. https://doi.org/10.3389/fcell.2020.00360

38. B. Stix, T. Kähne, K. Sletten, J. Raynes, A. Roessner, C. Röacken, Proteolysis of aa amyloid fibril proteins by matrix metalloproteinases -1, -2, and -3, *Am. J. Pathol.*, **159** (2001), 561–570. https://doi.org/10.1016/S0002-9440(10)61727-0
39. M. C. Liao, W. E. V. Nostrand, Degradation of soluble and fibrillar amyloid β-protein by matrix metalloproteinase (mt1-mmp) in vitro, *Biochemistry*, **49** (2010), 1127–1136.
https://doi.org/10.1021/bi901994d

40. D. Madsen, T. Bugge, The source of matrix-degrading enzymes in human cancer: problems of research reproducibility and possible solutions, *J. Cell Biol.*, **209** (2015), 195–198.
https://doi.org/10.1083/jcb.201501034

41. M. Aristorena, E. Gallardo-Vara, M. Vicen, M. de Las Casas-Engel, L. Ojeda-Fernandez, C. Nieto, et al., MMP-12, secreted by pro-inflammatory macrophages, targets endoglin in human macrophages and endothelial cells, *Int. J. Mol. Sci.*, **20** (2019), 3107.
https://doi.org/10.3390/ijms20123107

42. W. C. Huang, G. B. Sala-Newby, A. Susana, J. L. Johnson, A. C. Newby, Classical macrophage activation up-regulates several matrix metalloproteinases through mitogen activated protein kinases and nuclear factor-kb, *PLoS One*, **7** (2012), 1–14.
https://doi.org/10.1371/journal.pone.0042507

43. P. Vitorino, T. Meyer, Modular control of endothelial sheet migration, *Genes Dev.*, **22** (2008), 3268–3281.
https://doi.org/10.1101/gad.1725808

44. L. E. Tracy, R. A. Minasian, E. J. Caterson, Extracellular matrix and dermal fibroblast function in the healing wound, *Adv. Wound Care*, **5** (2016), 119–136.
https://doi.org/10.1089/wound.2014.0561

45. A. Alsisi, R. Eftimie, D. Trucu, Non-local multiscale approach for the impact of go or grow hypothesis on tumour-viruses interactions, *Math. Biosci. Eng.*, **18** (2021), 5252–5284.
https://doi.org/10.3934/mbe.2021267

46. A. Alsisi, R. Eftimie, D. Trucu, Non-local multiscale approaches for tumour-oncolytic viruses interactions, *Math. Appl. Sci. Eng.*, **1** (2020), 249–273.
https://doi.org/10.5206/mase/10773

47. Y. Koyama, K. Norose-Toyoda, S. Hirano, M. Kobayashi, M. Ebihara, I. Someki, et al., Type I collagen is a non-adhesive extracellular matrix for macrophages, *Arch. Histol. Cytol.*, **63** (2000), 71–79.
https://doi.org/10.1679/aohc.63.71

48. J. Y. Hsieh, M. T. Keating, T. D. Smith, V. S. Meli, E. L. Botvinick, W. F. Liu, Matrix crosslinking enhances macrophage adhesion, migration, and inflammatory activation, *APL Bioeng.*, **3** (2019), 016103.
https://doi.org/10.1063/1.5067301

49. G. F. Weber, M. A. Bjerke, D. W. DeSimone, Integrins and cadherins join forces to form adhesive networks, *J. Cell Sci.*, **124** (2011), 1183–1193.
https://doi.org/10.1242/jcs.064618

50. J. M. Teddy, J. M. Kulesa, In vivo evidence for short- and long-range cell communication in cranial neural crest cells, *Development (Cambridge, England)*, **131** (2004), 6141–6151.
https://doi.org/10.1242/dev.01534

51. A. Gerisch, A. Chaplain, Mathematical modelling of cancer cell invasion of tissue: local and non-local models and the effect of adhesion, *J. Theor. Biol.*, **250** (2008), 684–704.
https://doi.org/10.1016/j.jtbi.2007.10.026

52. T. J. R. Huges, *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*, Prentice Hall, Englewood Cliffs, New Jersey 07632, 1987.
53. P. Domschke, D. Trucu, A. Gerisch, M. Chaplain, Mathematical modelling of cancer invasion: implications of cell adhesion variability for tumour infiltrative growth patterns, *J. Theor. Biol.*, **361** (2014), 41–60. https://doi.org/10.1016/j.jtbi.2014.07.010

54. R. Shuttleworth, R. Trucu, Multiscale modelling of fibres dynamics and cell adhesion within moving boundary cancer invasion, *Bull. Math. Biol.*, **81** (2019), 2176–2219. https://doi.org/10.1007/s11538-019-00598-w

55. S. Luo, M. Benathan, W. Raffoul, R. G. Panizzon, D. V. Egloff, Abnormal balance between proliferation and apoptotic cell death in fibroblasts derived from keloid lesions, *Plast. Reconstr. Surg.*, **107** (2001), 87–96. https://doi.org/10.1023/A:1011941121102

56. S. Chhabra, N. Chhabra, A. Kaur, N. Gupta, Wound healing concepts in clinical practice of OMFS, *J. Maxillofac. Oral Surg.*, **250** (2017), 403–423. https://doi.org/10.1007/s12663-016-0880-z

57. G. Peyret, R. Mueller, J. d'Alessandro, J. Begnaud, P. Marcq, R. M. Mège, et al., Sustained oscillations of epithelial cell sheets, *Biophys. J.*, **117** (2019), 464–478. https://doi.org/10.1016/j.bpj.2019.06.013

58. R. B. Diller, A. J. Tabor, Role of the extracellular matrix (ecm) in wound healing: A review, *Biomimetics*, **7** (2022), 87. https://doi.org/10.3390/biomimetics7030087

59. J. Larouche, S. Sheoran, K. Maruyama, M. M. Martino, Immune regulation of skin wound healing: mechanisms and novel therapeutic targets, *Adv. Wound Care*, **7** (2018), 209–231. https://doi.org/10.1089/wound.2017.0761

60. J. Pang, M. Maienschein-Cline, M. Koh, Monocyte/macrophage heterogeneity during skin wound healing in mice, *J. Immunol.*, **209** (2022), 1999–2011. https://doi.org/10.4049/jimmunol.2200365

61. H. E. Talbott, S. Mascharak, M. Griffin, D. C. Wan, M. T. Longaker, Wound healing, fibroblast heterogeneity, and fibrosis, *Cell Stem Cell*, **29** (2022), 1161–1180. https://doi.org/10.1016/j.stem.2022.07.006

62. M. Alwuthaynani, R. Eftimie, D. Trucu, Inverse problem approaches for mutation laws in heterogeneous tumours with local and nonlocal dynamics, *Math. Biosci. Eng.*, **19** (2022), 3720–3747. https://doi.org/10.3934/mbe.2022171

63. J. van Kan, A. Segal, F. Vermolen, *Numerical Methods in Scientific Computing*, Delft Academic Press, Mekelweg 4 2628 CD Delft, Netherlands, 2014.

64. H. P. Langtangen, A FEniCS tutorial, in *Automated Solution of Differential Equations by the Finite Element Method, The FEniCS Book* (eds. A. Logg, K. A. Mardal, G. N. Wells), Springer, Berlin Heidelberg, (2012), 1–73. https://doi.org/10.1007/978-3-642-23099-8_1
# Appendix

## A. Reduction of the non-local terms

**Gaussian kernel:** Let $T_1, T_2$ represent the components of the term under the integral sign in Eq (2.16) namely:

$$T_i := \int_{-R}^{R} \int_{-R}^{R} \frac{y_i^2}{2\pi\sigma^2\sqrt{y_1^2 + y_2^2}} \exp\left(-\frac{y_1^2 + y_2^2}{2\sigma^2}\right) dy_1 dy_2, \quad i = 1, 2.$$

Thus, for $T_1$ we have

$$ 
\begin{aligned}
T_1 &= \int_{-R}^{R} \int_{-R}^{R} \frac{|y_1|}{2\pi\sigma^2} \cdot \frac{|y_2|}{\sqrt{|y_1|^2 + |y_2|^2}} \cdot \underbrace{\exp\left(-\frac{y_1^2 + y_2^2}{2\sigma^2}\right)}_{\le 1} dy_1 dy_2 \\
&\le \int_{-R}^{R} \int_{-R}^{R} \frac{|y_1|}{2\pi\sigma^2} dy_1 dy_2 = \int_{-R}^{R} \mathcal{L} \cdot \int_{0}^{R} \frac{y_1}{2\pi\sigma^2} dy_1 dy_2 \\
&= \int_{-R}^{R} \frac{R^2}{2\pi\sigma^2} dy_2 = \frac{R^3}{\pi\sigma^2}
\end{aligned}
\quad (4.1) $$

Following identical steps (as for $T_1$) and changing the order of integration, we obtain that the same upper bound holds true also for $T_2$.

Therefore $\lim_{R\to 0} T_{1,2}(R) = 0$, and the non-local fluxes disappear when we localise the integrals with Gaussian kernels.

**Cone-shaped kernel:** Substituting Eq (2.6) in Eq (2.20) gives the following local approximation of the non-local term:

$$ 
\begin{aligned}
A_{(\cdot)}[g, f, m, e] &= \frac{1}{R} \int_{0}^{R} r \int_{0}^{2\pi} \mathbf{n}(\theta) K(r) G(u(\mathbf{x}, t)) \, d\theta dr \\
&\quad + \frac{1}{R} \int_{0}^{R} r \int_{0}^{2\pi} \mathbf{n}(\theta) K(r) \langle \nabla_u G(u(\mathbf{x}, t)) \nabla u(\mathbf{x}, t), r\mathbf{n}(\theta) \rangle \, d\theta dr \\
&= \frac{1}{R} \int_{0}^{R} r \frac{3}{\pi R^2} \left(1 - \frac{r}{R}\right) dr dr \underbrace{\int_{0}^{2\pi} (\cos\theta, \sin\theta) \, d\theta}_{=(0,0)} \\
&\quad + \frac{\pi}{R} \nabla_u G(u(\mathbf{x}, t)) \nabla u(\mathbf{x}, t) \int_{0}^{R} r^2 \frac{3}{\pi R^2} \left(1 - \frac{r}{R}\right) dr \\
&= \frac{1}{4} \nabla_u G(u(\mathbf{x}, t)) \nabla u(\mathbf{x}, t)
\end{aligned}
$$
## B. Cell growth: logistic vs. truncated logistic

Because experimental studies [55] have shown that keloid fibroblasts have lower death rates and reach higher densities (and implicitly cell death due to overcrowding was reduced), we decided to investigate the impact of replacing a classical logistic cell growth with a truncated logistic growth (see Eq (3.10)). In Figure A1 we show that the solution of a simple logistic growth equation ($df/dt = p_f f(1 - \rho(u))$, $dm/dt = p_m m(1 - \rho(u))$) is the same as the solution of a truncated logistic growth ($df/dt = p_f f(1 - \rho(u))^+$, $dm/dt = p_m m(1 - \rho(u))^+$) if the initial condition is below the carrying capacity (as are usually our conditions). Note that an initial condition above the carrying capacity does not lead to a reduction in cell population size (since there is no cell death due to overcrowding). Therefore, in the situation where the spatial flux leads to local cell population overcrowding, the truncated logistic does not reduce the size of that population (while no further population is added).

<figure><img src="image_19.png" /><figcaption>**Figure A1.** Solution of ODE equations for temporal cell growth: $df/dt = p_f f(1 - \rho(u))$ (red continuous curve) and $df/dt = p_f f(1 - \rho(u))^+$ (blue dashed curve), for different initial conditions.</figcaption></figure>

## C. FEM discretization of non-local and local models

Using the backward (Implicit) Euler method as in Eq (3.5) for the time discretization, we have: find $\{c_{\tau}^g\}_{\tau=1}^l, \{c_{\tau}^f\}_{\tau=1}^l, \{c_{\tau}^m\}_{\tau=1}^l, \{c_{\tau}^e\}_{\tau=1}^l$ such that

$$ 
\begin{aligned}
\sum_{\tau=0}^{l} \frac{c_{\tau}^{g,N+1} - c_{\tau}^{g,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{j}(\mathbf{x}) = -D_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{j}(\mathbf{x}) \\
- \lambda_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \psi_{\tau}(\mathbf{x}) \psi_{j}(\mathbf{x}) + \left( p_g^f \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \right. \\
\left. + p_g^m \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \right) \psi_{j}(\mathbf{x}) \\
+ D_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \frac{\partial}{\partial n} \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \psi_{\tau}(\mathbf{x}) \psi_{j}(\mathbf{x}), \quad \forall j \in \{1, \dots, l\}, \quad (4.2a)
\end{aligned}
$$
$$
\begin{aligned}
\sum_{\tau=0}^{l} \frac{c_{\tau}^{f,N+1} - c_{\tau}^{f,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) = -D_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{o}(\mathbf{x}) + \nabla \cdot \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) A_f^N \right) \psi_o(\mathbf{x}) - \lambda_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) \\
+ p_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \in \psi_o(\mathbf{x}) \\
+ D_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \frac{\partial}{\partial n} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}), \quad \forall o \in \{1, \dots, l\},
\end{aligned}
\quad (4.2b)
$$

$$
\begin{aligned}
\sum_{\tau=0}^{l} \frac{c_{\tau}^{m,N+1} - c_{\tau}^{m,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \, d\mathbf{x} = -D_{m} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{h}(\mathbf{x}) \, d\mathbf{x} + \nabla \cdot \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) A_{m}^{N} \right) \psi_{h}(\mathbf{x}) \\
- \lambda_{m} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) + p_{m} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \in \psi_{h}(\mathbf{x}) \\
+ D_{m} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \frac{\partial}{\partial n} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}), \quad \forall h \in \{1, \dots, l\},
\end{aligned}
\quad (4.2c)
$$

$$
\begin{aligned}
\sum_{\tau=0}^{l} \frac{c_{\tau}^{e,N+1} - c_{\tau}^{e,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{\gamma}(\mathbf{x}) = - \sum_{\tau=0}^{l} c_{\tau}^{e,N+1} \nabla \psi_{\tau}(\mathbf{x}) \left( \alpha_{f} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + \alpha_{m} \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \right) \psi_{\gamma}(\mathbf{x}) \\
+ \sum_{\tau=0}^{l} c_{\tau}^{e,N+1} \psi_{\tau}(\mathbf{x}) p_{e} \in \psi_{\gamma}(\mathbf{x}), \quad \forall \gamma \in \{1, \dots, l\},
\end{aligned}
\quad (4.2d)
$$

Next, we formulate a barycentral approximation for the non-local terms $A_{f,m}^N(\mathbf{x})$. We show only Case II. (macrophage-ECM adhesion present: $S_{me} > 0$), since Case I. (no macrophage-ECM adhesion) is obtained for $S_{me} = 0$.

$$
A_f^N(\mathbf{x}) = \sum_{T \in \mathcal{A}(\mathbf{x})} K(\|\mathbf{x} - \mathbf{q}(T)\|_2) \sum_{\tau=0}^{l} \left( S_{ff,\tau}^N(\mathbf{q}(T)) c_{\tau}^{f,N} + \right. \\
\left. \qquad S_{fm,\tau}^N(\mathbf{q}(T)) c_{\tau}^{m,N} + S_{fe,\tau}^N(\mathbf{q}(T)) c_{\tau}^{e,N} \right) \psi_{\tau}(\mathbf{q}(T)) \left( 1 - \sum_{\tau=0}^{l} \left( c_{\tau}^{g,N} + c_{\tau}^{f,N} \right. \right. \\
\left. \left. + c_{\tau}^{m,N} + c_{\tau}^{e,N} \right) \psi_{\tau}(\mathbf{q}(T)) \right) \chi_T(\mathbf{q}(T)), \quad (4.3a)
$$

$$
A_m^N(\mathbf{x}) = \sum_{T \in \mathcal{A}(\mathbf{x})} K(\|\mathbf{x} - \mathbf{q}(T)\|_2) \sum_{\tau=0}^{l} \left( S_{mm,\tau}^N(\mathbf{q}(T)) c_{\tau}^{m,N} + \right. \\
\left. \qquad S_{mf,\tau}^N(\mathbf{q}(T)) c_{\tau}^{m,N} + S_{me,\tau}^N(\mathbf{q}(T)) c_{\tau}^{e,N} \right) \psi_{\tau}(\mathbf{q}(T)) \left( 1 - \sum_{\tau=0}^{l} \left( c_{\tau}^{g,N} + c_{\tau}^{f,N} + c_{\tau}^{m,N} \right. \right. \\
\left. \left. + c_{\tau}^{e,N} \right) \psi_{\tau}(\mathbf{q}(T)) \right) \chi_T(\mathbf{q}(T)), \quad (4.3b)
$$

Here $\mathcal{A}(\mathbf{x}) := \{T \in \mathcal{T} : \|\mathbf{x} - \mathbf{q}(T)\|_2 < R\}$ is the barycentral approximation, while $\mathbf{q} : \mathcal{T} \to \mathbb{R}$ gives the barycenter coordinates of triangles $T \in \mathcal{T}$ i.e.,

$$
\mathbf{q}(T) := \text{barycenter}(T), \quad \forall T \in \mathcal{T}
$$
### C1. FEM discretisation for local model (2.17) – reduced from the non-local model with Gaussian Kernel

The weak form of our coupled dynamics shown in Eq (2.17) can now be restated in terms of the basis functions [63] as follows:

find {$c_{\tau}^g$}$_{\tau=1}^l$, {$c_{\tau}^f$}$_{\tau=1}^l$, {$c_{\tau}^m$}$_{\tau=1}^l$, {$c_{\tau}^e$}$_{\tau=1}^l$ such that

$$ \sum_{\tau=0}^{l} \frac{c_{\tau}^{g,N+1} - c_{\tau}^{g,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{j}(\mathbf{x}) = -D_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{j}(\mathbf{x}) - \lambda_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \psi_{\tau}(\mathbf{x}) \psi_{j}(\mathbf{x}) + \left( p_g^f \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + p_g^m \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \right) \psi_{j}(\mathbf{x}) + D_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot n \psi_{j}(\mathbf{x}), \quad \forall j \in \{1, \dots, l\}, \quad (4.4a) $$

$$ \sum_{\tau=0}^{l} \frac{c_{\tau}^{f,N+1} - c_{\tau}^{f,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) = -D_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{o}(\mathbf{x}) \\ - \lambda_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) + p_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \epsilon \psi_{o}(\mathbf{x}) \\ + D_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \frac{\partial}{\partial n} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}), \quad \forall o \in \{1, \dots, l\}, \quad (4.4b) $$

$$ \sum_{\tau=0}^{l} \frac{c_{\tau}^{m,N+1} - c_{\tau}^{m,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) = -D_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{h}(\mathbf{x}) - \lambda_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\ + p_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \epsilon \psi_{h}(\mathbf{x}) \\ + D_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \frac{\partial}{\partial n} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}), \quad \forall h \in \{1, \dots, l\}, \quad (4.4c) $$

$$ \sum_{\tau=0}^{l} \frac{c_{\tau}^{e,N+1} - c_{\tau}^{e,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{y}(\mathbf{x}) = - \sum_{\tau=0}^{l} c_{\tau}^{e,N+1} \nabla \psi_{\tau}(\mathbf{x}) \left( \alpha_f \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + \alpha_m \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \right) \psi_{y}(\mathbf{x}) \\ + \sum_{\tau=0}^{l} c_{\tau}^{e,N+1} \psi_{\tau}(\mathbf{x}) p_e \epsilon \psi_{y}(\mathbf{x}), \quad \forall \gamma \in \{1, \dots, l\}. \quad (4.4d) $$

Here, $S_{\kappa} = S_{\kappa}^{\max} \vartheta$, $\kappa \in \{ff, fm, mf, mm, fe, me\}$ represents the strength of cell-cell and cell-ECM interactions discretized in time and space.

$$ \nabla S_{\kappa} = S_{\kappa}^{\max} (\varepsilon - \vartheta^2) \left( \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right), \quad \kappa \in \{ff, fm, mf, mm, fe, me\}, \quad (4.5) $$

represents their gradients, where
$$ \epsilon = \frac{1}{1 + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x})}, \quad (4.6) $$

$$ \vartheta = \frac{\sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x})}{1 + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x})}, \quad (4.7) $$

and $S_g^\kappa = \frac{\partial}{\partial g} S_\kappa = S_e^\kappa = \frac{\partial}{\partial e} S_\kappa = \frac{1}{e^2}$, $\kappa \in \{ff, fm, mf, mm, fe, me\} \equiv \{1, 2, 3, 4, 5, 6\}$. Also, we denote by $\epsilon$ the overcrowding term $1 - \rho(u)$ discretized in time and space, i.e.,

$$ \epsilon = 1 - \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) - \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) - \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) - \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \quad (4.8a) $$

## C2. FEM discretisation for local model (2.22) – reduced from non-local model with cone-shaped kernel

Below, we show the FEM discretisation for the local models. However, since Case I (no macrophage-ECM adhesion: $S_{me} = 0$) is just a simplification of Case II (that includes macrophage-ECM adhesion), in the following, we present only the FEM discretisation for Case II ($S_{me} > 0$).

The weak form of system (2.22) (obtained via Taylor-series expansion of the non-local model with the cone-shaped kernel (2.6)) is given as follows in terms of basis functions:

Find $\{c_{\tau}^g\}_{\tau=1}^l, \{c_{\tau}^f\}_{\tau=1}^l, \{c_{\tau}^m\}_{\tau=1}^l, \{c_{\tau}^e\}_{\tau=1}^l$ such that

$$ \sum_{\tau=0}^{l} \frac{c_{\tau}^{g,N+1} - c_{\tau}^{g,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{j}(\mathbf{x}) = -D_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{j}(\mathbf{x}) - \lambda_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \psi_{\tau}(\mathbf{x}) \psi_{j}(\mathbf{x}) + \left( p_g^f \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + p_g^m \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \right) \psi_{j}(\mathbf{x}) + D_g \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot n \psi_{j}(\mathbf{x}), \quad \forall j \in \{1, \dots, l\}, \quad (4.9a) $$

$$ \sum_{\tau=0}^{l} \frac{c_{\tau}^{f,N+1} - c_{\tau}^{f,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) = D_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{o}(\mathbf{x}) + \frac{1}{8} \mu_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot A_f^N(\mathbf{x}) \\ + \mu_f \frac{1}{4} \left( (\nabla S_{ff}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) \right) \\ - w_g \left( \nabla S_{ff} \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) $$
$$
\begin{aligned}
& + \mathbf{S}_{ff} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) \\
& - 2w_f \Bigg( \nabla \mathbf{S}_{ff} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{ff} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) \\
& - w_m \Bigg( \nabla \mathbf{S}_{ff} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{ff} \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) \\
& - w_e \Bigg( \nabla \mathbf{S}_{ff} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{ff} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) \\
& - w_f \Bigg( \mathbf{S}_{fm} \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \nabla \mathbf{S}_{fm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla \mathbf{S}_{fe} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{fe} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) \\
& + S_g^1 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^1 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_{o}(\mathbf{x}) \\
& - w_g \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^1 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) S_g^1 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \Bigg)^2 S_g^1 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \psi_{o}(\mathbf{x}) \\
& - w_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \nabla S_g^1 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg)
\end{aligned}
$$
$$
\begin{aligned}
& + 2S_g^1 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x})\psi_o(\mathbf{x}) \\
& - w_m \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^1 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) S_g^1 \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) S_g^1 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x})\psi_o(\mathbf{x})dx \\
& - w_e \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla\psi_{\tau}(\mathbf{x}) \nabla S_g^1 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) S_g^1 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) S_g^1 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x})\psi_o(\mathbf{x})dx \\
& - w_g \Bigg( S_{ff} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_{ff} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) + S_{fm} \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) S_{fm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x})\psi_o(\mathbf{x}) + \\
& \qquad \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla\psi_{\tau}(\mathbf{x}) \nabla S_g^2 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \\
& \qquad + S_g^2 \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x})\psi_o(\mathbf{x}) \\
& - w_g \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_a(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^2 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) S_g^2 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) S_g^2 \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x})\psi_o(\mathbf{x})dx \\
& - w_f \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^2 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla\psi_{\tau}(\mathbf{x}) \Bigg)
\end{aligned}
$$
$$
\begin{aligned}
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) S_g^2 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) + \\
& \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) S_g^2 \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) \right. \\
& - w_m \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^2 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right. \right. \\
& \left. + 2S_g^2 \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) \right) \\
& - w_e \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^2 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \left. + \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) S_g^2 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \left. \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) S_g^2 \sum_{\tau=0}^{l} c_{\tau}^{m,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) \right) \\
& + \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \nabla S_g^5 \cdot \nabla \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \right. \\
& + S_g^5 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \left) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) \\
& - w_g \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}), \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^5 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \left. + S_g^5 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \left. S_g^5 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) \right) \\
& - w_f \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^5 \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \left. + S_g^5 \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& + S_g^5 \sum_{\tau=0}^{l} c_{\tau}^{e,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) S_e^2 \epsilon + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) S_e^5 \\
& + \epsilon + S_{fe} \epsilon - w_e S_{fe} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_o(\mathbf{x}) \right)
\end{aligned}
$$
$$
\begin{aligned}
&+ \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) S_e^2 \epsilon - w_e S_{ff} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \right. \\
&\quad - w_e S_{fm} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) S_e^2 \epsilon + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) S_e^5 \epsilon + S_{fe} \epsilon \\
&\quad \left. - w_e S_{fe} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right) \frac{\partial}{\partial n} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) \\
&\quad - \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \left( S_{fm} \epsilon - w_m S_{fm} \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) - w_m S_{ff} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \right. \\
&\quad \left. - w_m S_{fe} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_o(\mathbf{x}) \\
&\quad + \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \left( S_{fm} \epsilon - w_m S_{fm} \sum_{\tau=0}^{l} c_{\tau}^{m,N} \psi_{\tau}(\mathbf{x}) - w_m S_{ff} \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \right. \\
&\quad \left. - w_m S_{fe} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \frac{\partial}{\partial n} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) \\
&\quad - \lambda_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \epsilon \psi_o(\mathbf{x}) \\
&\quad + D_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \frac{\partial}{\partial n} \sum_{\tau=0}^{l} c_{\tau}^{f} \psi_{\tau}(\mathbf{x}) \psi_o(\mathbf{x}), \quad \forall o \in \{1, \dots, l\}, \quad (4.9b)
\end{aligned}
$$

$$
\begin{aligned}
\sum_{\tau=0}^{l} \frac{c_{\tau}^{m,N+1} - c_{\tau}^{m,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{\tau}(\mathbf{x}) = & -D_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{\tau}(\mathbf{x}) \\
& + \frac{1}{8} \mu_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \left( \left( S_{mf} \epsilon - w_f S_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \right. \right. \\
& \left. \left. - w_f S_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) - w_f S_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \left. + \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \epsilon S_g^4 - w_g S_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) - w_g S_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \right. \right. \\
& \left. \left. + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \epsilon S_g^3 + S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \epsilon - w_g S_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \left. + \left( S_{mm} \epsilon - w_m S_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) - w_m S_{ff} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \right. \right. \\
& \left. \left. - w_m S_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) + \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \epsilon S_e^4 \right)
\end{aligned}
$$
$$
\begin{aligned}
& -w_e \mathbf{S}_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) - w_e \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \in S_e^3 \\
& + S_e^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) - w_e \mathbf{S}_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \Biggr) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla\psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& + \mu_m \frac{1}{2} \Biggl( \nabla \mathbf{S}_{mf} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \Biggr) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_g \Biggl( \nabla \mathbf{S}_{mf} \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \Biggr) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - 2w_f \Biggl( \nabla \mathbf{S}_{mf} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \Biggr) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_m \Biggl( \nabla \mathbf{S}_{mf} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \Biggr) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_e \Biggl( \nabla \mathbf{S}_{mf} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \Biggr) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_f \mathbf{S}_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \Biggl( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla \mathbf{S}_{mm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla \mathbf{S}_{me} \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \\
& + \mathbf{S}_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla\psi_{\tau}(\mathbf{x}) \Biggr) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& + S_g^4 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla\psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla\psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla\psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x})
\end{aligned}
$$
$$
\begin{aligned}
& -w_g \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_g^4 \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) S_g^4 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& \qquad + 2S_g^4 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_f \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_g^4 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_g^4 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_e \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_g^4 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) S_g^4 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_g \left( S_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_{mm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) + S_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_{mf} \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& + \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^3 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& + S_g^3 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x})
\end{aligned}
$$
$$
\begin{aligned}
& -w_g \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^3 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_g^3 \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) S_g^3 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& \qquad - w_f \left( \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \right)^2 \nabla S_g^3 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. 2 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_g^3 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& \qquad - w_m \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^3 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_g^3 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_g^3 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& \qquad - w_e \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^3 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) S_g^3 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_g^3 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^6 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& \qquad - w_g \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^6 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \\
& \qquad + S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x})
\end{aligned}
$$
$$
\begin{aligned}
& -w_f \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_g^6 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_m \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_g^6 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_e \left( \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right)^2 \nabla S_g^6 \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& 2S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_g \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_{me} \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + S_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& + \nabla S_{mm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_g \left( S_{mm} \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_{mm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_f \left( S_{mm} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_{mm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) \\
& - w_e \left( S_{mm} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \right)
\end{aligned}
$$
$$
\begin{aligned}
& + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla \mathbf{S}_{mm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) - \\
& \qquad 2w_m \Bigg( \mathbf{S}_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla \mathbf{S}_{mm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_m \Bigg( \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla \mathbf{S}_{mf} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_m \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla \mathbf{S}_{me} \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + \mathbf{S}_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad + \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_e^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + S_e^4 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_g \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_e^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_e^4 \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) S_e^4 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_f \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_e^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_e^4 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_e^4 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_m \Bigg( \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \Bigg)^2 \nabla S_e^4 \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) +
\end{aligned}
$$
$$
\begin{aligned}
& 2 \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_{e}^{4} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \right) \mathrm{d}\mathbf{x} \\
& - w_{e} \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{4} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_{e}^{4} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) S_{e}^{4} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \right) \\
& - w_{e} \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_{mm} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + S_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \right) \\
& - w_{e} \left( S_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_{mf} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \right) \\
& + \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{3} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + S_{e}^{3} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \right) \\
& - w_{g} \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{3} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right. \\
& \qquad \left. + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_{e}^{3} \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) S_{e}^{3} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \right) \\
& - w_{f} \left( \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \right)^{2} \nabla S_{e}^{3} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. 2 \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_{e}^{3} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \right) \\
& - w_{m} \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{3} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right)
\end{aligned}
$$
$$
\begin{aligned}
&+ \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_{e}^{3} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \\
&\sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_{e}^{3} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
&- w_{e} \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{3} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
&+ \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_{e}^{3} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg) \\
&+ \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) S_{e}^{3} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
&+ \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{6} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
&+ S_{e}^{6} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
&- w_{g} \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{6} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
&+ S_{e}^{6} \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \\
&S_{e}^{6} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
&- w_{f} \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{6} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
&+ S_{e}^{6} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg) \\
&+ S_{e}^{6} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
&- w_{m} \Bigg( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_{e}^{6} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \\
&+ S_{e}^{6} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \Bigg) \\
&+ S_{e}^{6} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \Bigg|_{\tau=0} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x})
\end{aligned}
$$
$$
\begin{aligned}
& -w_e \left( \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right)^2 \nabla S_e^6 \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. 2S_e^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad + \left( \nabla S_{me} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_g \left( \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x}) \nabla S_{me} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. + S_{me} \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_f \left( \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \nabla S_{me} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. + S_{me} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_m \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \nabla S_{me} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. + S_{me} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - w_e \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \nabla S_{me} \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) + \right. \\
& \qquad \left. + S_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_{h}(\mathbf{x}) \\
& \qquad - \mu_f \frac{1}{2} \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( S_{mm} \epsilon - w_f \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \right) + \right. \\
& \qquad \left. - w_f \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_{mm} - w_f S_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right) \sum_{\tau=0}^{l} c_{\tau}^{f,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{h}(\mathbf{x}) \\
& \qquad + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( S_{mf} \epsilon - w_f \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) - w_f \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_{mm} - w_f S_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \right) \frac{\partial \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x})}{\partial n} \psi_{h}(\mathbf{x}) \\
& \qquad - \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_g^4 \epsilon - w_g S_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \right)
\end{aligned}
$$
$$
\begin{aligned}
& - w_g \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_g^3 \epsilon + S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \epsilon \\
& - w_g \mathbf{S}_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{g,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{h}(\mathbf{x}) \right) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_g^4 \epsilon - w_g \mathbf{S}_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \right) \\
& - w_g \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_g^3 \epsilon + S_g^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \epsilon \\
& - w_g \mathbf{S}_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \left( \frac{\partial \sum_{\tau=0}^{l} c_{\tau}^{g,N} \psi_{\tau}(\mathbf{x})}{\partial n} \right) \psi_h(\mathbf{x}) \\
& - \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_e^4 \epsilon - w_e \mathbf{S}_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \right) \\
& - w_e \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_e^3 \epsilon + S_e^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \epsilon + \mathbf{S}_{me} \epsilon \\
& - w_e \mathbf{S}_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{e,N} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{h}(\mathbf{x}) \right) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) S_e^4 \epsilon - w_e \mathbf{S}_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \right) \\
& - w_e \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) S_e^3 \epsilon + S_e^6 \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \epsilon + \mathbf{S}_{me} \epsilon \\
& - w_e \mathbf{S}_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \left( \frac{\partial \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x})}{\partial n} \right) \psi_h(\mathbf{x}) \\
& - \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( \mathbf{S}_{mm} \epsilon w_m \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) - w_m \mathbf{S}_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \right) \\
& - w_m \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \left( \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \nabla \psi_{\tau}(\mathbf{x}) \cdot \nabla \psi_{h}(\mathbf{x}) \right) \\
& + \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \left( \mathbf{S}_{mm} \epsilon - w_m \mathbf{S}_{mf} \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) - w_m \mathbf{S}_{mm} \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \right) \\
& - w_m \mathbf{S}_{me} \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) \left( \frac{\partial \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x})}{\partial n} \right) \psi_h(\mathbf{x}) \\
& - \lambda_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \psi_h(\mathbf{x}) + p_m \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x}) \epsilon \psi_h(\mathbf{x})
\end{aligned}
$$
$$ + D_m \frac{\partial \sum_{\tau=0}^{l} c_{\tau}^{m,N+1} \psi_{\tau}(\mathbf{x})}{\partial n} \psi_h(\mathbf{x}), \quad (4.9c) $$

$$ \sum_{\tau=0}^{l} \frac{c_{\tau}^{e,N+1} - c_{\tau}^{e,N}}{\Delta t} \psi_{\tau}(\mathbf{x}) \psi_{\gamma}(\mathbf{x}) = -\sum_{\tau=0}^{l} c_{\tau}^{e,N+1} \nabla \psi_{\tau}(\mathbf{x}) \left( \alpha_f \sum_{\tau=0}^{l} c_{\tau}^{f,N+1} \psi_{\tau}(\mathbf{x}) \right. \\ \left. + \alpha_m \sum_{\tau=0}^{l} c_{\tau}^{f,N} \psi_{\tau}(\mathbf{x}) \right) \psi_{\gamma}(\mathbf{x}) + \sum_{\tau=0}^{l} c_{\tau}^{e,N} \psi_{\tau}(\mathbf{x}) p_e \in \psi_{\gamma}(\mathbf{x}), \quad \forall \gamma \in \{1, \dots, l\}. \quad (4.9d) $$

### D. Time-marching in FEniCS

We re-write Eq (3.5) in the standard form $F(u^{N+1}; u^N; v) = 0$, where

$$ F(u^{N+1}; u^N; v) = \int_{\Omega} (u^{N+1} - u^N - \Delta t G(u^{N+1})) v \, dx. \quad (4.10) $$

Substituting Eq (3.4) into Eq (4.10) we have:

$$ F(c_{\tau}^{u_i^{N+1}}; c_{\tau}^{u_i^N}; \psi_{\tau}(\mathbf{x})) = \sum_{\tau=0}^{l} (c_{\tau}^{u_i^{N+1}} - c_{\tau}^{u_i^N}) \psi_{\tau}(\mathbf{x}) - \sum_{\tau=0}^{l} \Delta t G(c_{\tau}^{u_i^{N+1}}; \psi_{\tau}(\mathbf{x})) \psi_{j}(\mathbf{x}) \, dx, \forall j \in \{1, \dots, l\} \quad (4.11) $$

Below is an algorithmic form of the FEniCS procedure at each time step for nonlinear F [64]:

set some stopping time $T$

$$ t = \Delta t $$

while $t \le T$ do

evaluate $J$ (i.e, the Jacobian matrix of $F$)

solve $F = 0$ using a nonlinear solver that relies on the Jacobian of $F$ (i.e., Newton's method)

$$ t \leftarrow t + \Delta t \quad (\text{update time}) $$

$$ z_{\text{prev}} \leftarrow z \quad (\text{update solution}). $$

end while

### E. Mesh Refinement Analysis: A justification for the choice of mesh size

In deciding the choice of mesh dimension for our simulations, a mesh refinement analysis was carried out using the fine mesh 256 × 256 as our exact solution while mesh sizes 4 × 4, 8 × 8, 16 × 16, 32 × 32, 64 × 64, and 128 × 128 represented our approximate solutions. Figure A2 is a plot of the L₂ errors obtained against the different mesh sizes. The errors observed using mesh sizes 32 × 32 are insignificant compared to those of mesh sizes 4 × 4 and 8 × 8. For this reason, we decided to use a mesh size of 32 × 32 for all simulations presented in this study.

### F. Simulation results for two new initial conditions

Next we show the simulation results obtained for non-local model Case II with cone-shaped with two new initial conditions: one describing a linear but irregular cut in the tissue, as given by Eq (4.12d),

$$ g(\mathbf{x}, 0) = 0.1, \quad (4.12a) $$
<figure><img src="image_20.png" /><figcaption>**Figure A2.** L₂ error of the fibroblast solution against mesh grid sizes {2ᵢ × 2ᵢ : i = 2, ..., 5}, calculated with respect to the solution at 128 × 128.</figcaption></figure>

$$f(\mathbf{x}, 0) = 0.4[(0.5 + 0.5 \tanh(20x_1 - 9)) + (0.5 + 0.5 \tanh(-20x_1 - 9))], \quad (4.12b)$$

$$m(\mathbf{x}, 0) = 0.1[(0.5 + 0.5 \tanh(20x_1 - 9)) + (0.5 + 0.5 \tanh(-20x_1 - 9))], \quad (4.12c)$$

$$e(\mathbf{x}, 0) = 1.0[(0.5 + 0.5 \tanh(20x_1 - 9)) + 0.5 + 0.5 \tanh(-20x_1 - 9) \\ + 0.5 - ((0.2 + 0.2 \tanh(20x_1 - 3)) - 0.2) - 0.2 \tanh(-20(x_1 - 0.2) - 2) \\ - 0.1 - (0.2 + 0.2 \tanh(20x_1 - 0.03)) - 0.1) + (0.2 - 0.2 \tanh(-20x_1 - 0.001)) - 3]. \quad (4.12d)$$

and the other describing a circular wound, as given by Eq (4.13d),

$$g(\mathbf{x}, 0) = 0.1, \quad (4.13a)$$

$$f(\mathbf{x}, 0) = 0.4 \left( 1 - \exp \left( - \frac{x_1^2 + x_2^2}{0.04} \right) \right), \quad (4.13b)$$

$$m(\mathbf{x}, 0) = 0.1 \left( 1 - \exp \left( - \frac{x_1^2 + x_2^2}{0.04} \right) \right), \quad (4.13c)$$

$$e(\mathbf{x}, 0) = 1 - \exp\left(-\frac{x_1^2 + x_2^2}{0.04}\right). \quad (4.13d)$$

In Figures A3 and A4, we see that at time t = 100 both of the wounds heal normally as the ECM approaches its maximum level while the fibroblast and macrophage tissue level approach their baseline level. However, unlike the healing in the previous figures (for the original initial conditions), here we see that just before the final healing, the ECM is increased in the wound region (above the density of the ECM in the surrounding tissue) and then decreases towards the density of the surrounding tissue. We note that such a transient behaviour could have been possible also for the simulations performed with the original initial conditions. However, since we did not show the solutions at every single time step, such a detailed behaviour was lost before and was randomly revealed through Figures A3 and A4. This suggest that a more thorough analytical and numerical investigation of transient model dynamics is necessary to understand the behaviour of these non-local mathematical models.
<figure><img src="image_21.png" /><figcaption>**Figure A3.** *Initial Conditions: circular wound (see Eq (4.13d)).* Simulations of the *non-local model* (2.1) with *classical logistic term for both cells and ECM*. We consider case II with a cone-shaped kernel. The rows correspond to the spatial distribution of growth factor (*g*), fibroblast (*f*), macrophages (*m*), and ECM (*e*) at time points *t* = 0, *t* = 2, *t* = 20, and *t* = 100. The parameter values are listed in Table 1.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure A4.** *Initial Conditions: irregular linear cut (see Eq (4.12d)).* Simulations of the *non-local model* (2.1) with *classical logistic term for both cells and ECM*. We consider case II with a cone-shaped kernel. The rows correspond to the spatial distribution of growth factor (*g*), fibroblast (*f*), macrophages (*m*), and ECM (*e*), at time points *t* = 0, *t* = 2, *t* = 20, and *t* = 100. The parameter values are listed in Table 1.</figcaption></figure>
G. Sensitivity test for the sensing radius

Since there are some small differences between the non-local models and the corresponding local
models (obtained in the limit R → 0), in the following we perform a numerical sensitivity test for
the shape of the solutions as we vary the cells sensing radius R. In Figure A5 we show a series of
space-slice of these solutions for R = 0.08, R = 0.1, and R = 0.13, at times t = 2, 20, 40, and 100. We
see that increasing the sensing radius, leads to a faster remodeling of the ECM and a faster decay of
fibroblasts, macrophages, and the growth factor.

<figure><img src="image_23.png" /><figcaption>**Figure A5.** Sensitivity of the non-local model (2.1) with *classical logistic term for both cells* and *ECM* to the cells’ sensing radius *R*. We consider case II with a cone-shaped kernel. The columns correspond to the spatial distribution (along *x* axis, and at *y* = 0) for the growth factor (*g*), fibroblast (*f*), macrophages (*m*), and ECM (*e*), at time points *t* = 2, *t* = 20, *t* = 40, and *t* = 100. The parameter values are listed in Table 1.</figcaption></figure>

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

AIMS Press

©2023 the Author(s), licensee AIMS Press. This
is an open access article distributed under the
terms of the Creative Commons Attribution License
(http://creativecommons.org/licenses/by/4.0)