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

1.   Introduction

Wound healing is a complex process through which an organism tries to restore the integrity of biological tissue following 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: (ⅰ) by primary intention, where the incision is narrow and the wound heals as its edges are brought close together; (ⅱ) by secondary intention, where the wound heals from the bottom of the wound up. A normal wound healing process (i.e., either by primary or secondary intention) follows a series of interconnected phases (i.e., hemostasis, inflammation, proliferation and remodelling) ^[3], which eventually returns the wounded tissue close to its original state. However, when there is 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; they thus remain an uncovered clinical need. Deciphering the overlapping phases involved in normal and abnormal wound healing is key to understanding how these types of fibrosis are triggered and evolve over time. In what follows, we briefly detail the 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 ECM 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,8,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 clean the wound and 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 protein deposition and connective tissue contraction through the differentiation of fibroblasts into myofibroblasts, and, finally wound re-epithelialisation. During the proliferation phase, fibroblasts are chemo-attracted and migrate (particularly in response to TGF-$\beta$); they then 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 Ⅲ collagen being replaced with newly secreted type Ⅰ collagen. Such remodelling is the result of a "synthesis-degradation" balance under the control of matrix metalloproteinases that cause collagen breakdown. These proteinases, and their inhibitors are secreted by fibroblasts and macrophages as well. TGF-$\beta$ can inhibit the secretion of these metalloproteinases leading to the accumulation of collagen fibres ^[7], as observed in both hypertrophic scars and keloids.

Since the complexity of the above mentioned 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,14,15,16]}, while others focus on abnormal wounds ^[17,18,19]. Models usually describe healing by primary intention ^[20], focus mainly on the biochemical interactions between cells (or between cells and ECM) and do not include the biomechanical forces involved in wound contraction. In contrast, models that focus on mechanical interactions ^[2,9,21,22] usually describe healing by secondary intention ^[20]. However, to the best of our knowledge, there are not many models that concurrently investigate the biological mechanisms behind healing by primary or secondary intentions in the context of understanding normal versus 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 was to shed new light on the biological mechanisms that could generate different types of wounds (i.e., normal or abnormal) that heal by primary or secondary intention. To this end, we have considered a theoretical (i.e., modelling and numerical) approach to investigate some of the inflammatory and biomechanical aspects involved in the different phases of wound healing discussed above. Here we focus on the roles of TGF-$\beta$, macrophages and fibroblasts during the inflammation, proliferation and remodelling steps. The non-local model developed here accounts for the non-local impact of biomechanical interactions between cells, as well as among cells and the ECM as recent experimental studies have shown long-distance interactions between cells ^[23,24,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 primary or secondary.

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-$\beta$, $g({\bf x}, t)$, the density of fibroblasts, $f({\bf x}, t)$, the density of macrophages, $m({\bf x}, t)$, and the density of the ECM, $e({\bf 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 can be described by the following equations:

These equations incorporate the following biological assumptions:

● The changes in the concentration of the growth factor (Eq (2.1a)) at a position $({\bf 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,29,30], as described by the term $W(f, m)$ below:

where $p^{^f}_{_g}$ is the production rate for the growth factor due to fibroblasts ^[28,29,30], and $p^{^m}_{_g}$ is the production rate for the growth factor due to macrophages ^[28].

● The changes in the density of fibroblasts (Eq (2.1b)) at a position $({\bf 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 the coefficient $\mu_{_f}$, which is a consequence of cell-cell and cell-ECM adhesion ^[32,33]). In addition, the change in fibroblast densities is the result of fibroblast apoptosis ^[28] (at a rate $\lambda_{_{f}}$) and fibroblast 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 from other cells ^[23,24,25], in Eq (2.1b), $A_{_{f}}\left[g, f, m, e\right]$ 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 as well as the role of growth factor in these processes. A detailed description of this non-local flux term is given below. Moreover,

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 ECM respectively. Therefore, the proliferation of cells is described by the logistic term $f(1-\rho(u))$ that models cell proliferation according to nutrient availability, which is consistent with some experimental studies ^[35].

● The changes in the density of macrophages (Eq (2.1c)) at a position $({\bf x}, t)$ are a result of the flux of macrophages ^[36] (which consists of a linear diffusion term with the coefficient $D_{_{m}}$ and a cell migration term with the coefficient $\mu_{_m}$, which is a consequence of cell-cell and cell-ECM adhesion ^[32,33]), macrophage 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}}\left[g, f, m, e\right]$ denotes a non-local spatial flux term describing the adhesion processes between the fibroblasts, macrophages and ECM as well as the role of TGF-$\beta$ in these adhesive interactions. The proliferation of macrophages is described by the logistic term $m(1-\rho(u))$.

● The changes in the density of the 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)$ which is assumed to depend linearly on the density of fibroblasts (i.e., since dermal fibroblasts are closely linked to the ECM as both producer and resident cells ^[44]). Finally, as in ^[45,46], we consider a logistic-type remodelling of the ECM.

The adhesive cell-cell and cell-ECM interactions between the cells distributed at ${\bf x}$ and the surrounding cells and ECM perceived over a ball-shaped sensing region ${\bf B}({\bf x}, R): = {\bf x}+{\bf B}(0, R)$ of radius $R > 0$ are expressed via the non-local terms:

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

where $\|\cdot\|_{_{2}}$ is the usual Euclidean norm. The kernel $\mathrm{K}(\cdot):\left[0, R\right]\xrightarrow{}\left[0, 1\right]$ is a radially symmetric kernel that gives the interaction range of cells (i.e., interactions between the reference cell at position ${\bf x}$ and the neighbours at ${\bf x+y}$). Examples of such kernels are as follows:

${\rm{a}}.$ Gaussian kernel (see Figure 1(a))

${\rm{b}}.$ Cone-shaped kernel (see Figure 1(b))

Moreover, the term $(1-\rho (u))^{^{+}} : = \text{max} \{(1-\rho (u)), 0\}$ is included to avoid overcrowding within the non-local interactions. Finally, the function ${\bf \Gamma}_{_{f, m}} ({\bf x}+{\bf y}, t)$ describes the type and magnitude of cell-cell and cell-ECM adhesive interactions between cells at position ${\bf x}$ and the neighbours at position ${\bf x+y}$. To define these ${\bf \Gamma}_{_{f, m}}$ {functions}, we assume that the fibroblasts are cocultured with macrophages on the ECM. Note here that some experimental studies observed that macrophages cannot adhere to the ECM (at least to that of type Ⅰ 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:

${\rm{I}}.$ No macrophage-ECM adhesion.

${\rm{II}}.$ Including macrophage-ECM adhesion ($S_{_{me}} > 0$ in ${\Gamma}_{_{m}}$).

In the above equations, we consider the strengths of fibroblast-fibroblast interactions $({\bf S}_{_{ff}})$, fibroblast-macrophage interactions (${\bf S}_{_{fm}}$), fibroblast-ECM interactions (${\bf S}_{_{fe}}$), macrophage-macrophage interactions (${\bf S}_{_{mm}}$), macrophage-fibroblast interactions (${\bf S}_{_{mf}}$) and macrophage-ECM interactions (${\bf S}_{_{me}}$). Since these interaction strengths depend on TGF-$\beta$ and the presence of the ECM ^[49], to define them, we consider a monotonically increasing Hill-type function that depends on "$e+g$" and satisfies ${\bf S}_{_{j}}(e, g)$ = 0 for $e = 0$ and $g = 0$:

The above function is used to describe cell movement against the gradient of the ECM and neighbouring cells in the presence of the growth factor (and it 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 these dynamics, in Section 3.3, we use the following function for cell-cell and cell-ECM adhesion:

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 these 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 cells ^[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$; 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

where $A_{_{(.)}}[g, f, m, e]\in \{A_{_{m}}[g, f, m, e], A_{_{f}}[g, f, m, e]\}$ and ${\bf y} = (y_{_{1}}, y_{_{2}})$. The function $G(u({\bf x} + {\bf y}, t))$ is defined by the following expressions:

● for $A_{_{f}}$, we have

● for $A_{_{m}}$, we have the following general expression:

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

The truncated local Taylor expansion of $G(u({\bf x+y}, t))$ around ${\bf x}: = (x_{_{1}}, x_{_{2}})$ is given as

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

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

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

We emphasise 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 the 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 did 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 following form:

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

● For the fibroblasts non-local flux term $A_{_{f}}\left[g, f, m, e\right]$, the cell-cell and cell-ECM interactions are given by

● For the macrophages non-local flux term $A_{_{m}}\left[g, f, m, e\right],$ the cell-cell and cell-ECM interactions are given by

As before, in the above equation, condition ${\bf S}_{_{me}} = 0$ corresponds to Case Ⅰ (i.e., no macrophage-ECM adhesion), while condition ${\bf S}_{_{me}} > 0$ corresponds to Case Ⅱ (i.e., macrophage-ECM adhesion present).

The Taylor series expansion of $G(u({\bf x}+r{\bf n}(\theta), t)) = (G\circ u)({\bf x}+r\textbf{ n}(\theta), t)$ around $r = 0$ is given as

where

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

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

In the limit $R\to 0$, the non-local flux term $A_{_{f, m}}\left[g, f, m, e\right]$ approach the local term $A^{^{0}}_{_{f, m}}(g, f, m, e)$, where

In Eq (2.21a), the term ${\bf S}_{_{ff}}\nabla f + {\bf S}_{_{fm}}\nabla m + {\bf S}_{_{fe}}\nabla e$ on the right-hand side describes fluxes due to fibroblast-fibroblast, fibroblast-macrophage and fibroblast-ECM adhesions, with strengths ${\bf S}_{_{ff}}, {\bf S}_{_{fm}}$ and ${\bf S}_{_{fe}}$, respectively. The term $f\left(\frac{\partial}{\partial g}{\bf S}_{_{ff}}\nabla g + \frac{\partial}{\partial e}{\bf S}_{_{ff}}\nabla e\right) + m\left(\frac{\partial}{\partial g}{\bf S}_{_{fm}}\nabla g + \frac{\partial}{\partial e}{\bf S}_{_{fm}}\nabla e\right) + e\left(\frac{\partial}{\partial g}{\bf S}_{_{fe}}\nabla g + \frac{\partial}{\partial e}{\bf S}_{_{fe}}\nabla e\right)$ describes the changes in the strengths of fibroblast-fibroblast, fibroblast-macrophage, and fibroblast-ECM adhesion with respect to the growth factor and ECM, along the positive gradients of the growth factor and ECM, respectively. Finally, the term $- ({\bf S}_{_{ff}}f + {\bf S}_{_{fm}}m + {\bf 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 ${\bf S}_{_{mf}}\nabla f + {\bf S}_{_{mm}}\nabla m$ on the right-hand side represents macrophage-fibroblast and macrophage-macrophage adhesion with strengths ${\bf S}_{_{mf}}$ and ${\bf S}_{_{mm}}$, respectively. The term $f\left(\frac{\partial}{\partial g}{\bf S}_{_{mf}}\nabla g + \frac{\partial}{\partial e}{\bf S}_{_{mf}}\nabla e\right) + m\left(\frac{\partial}{\partial g}{\bf S}_{_{mm}}\nabla g + \frac{\partial}{\partial e}{\bf S}_{_{mm}}\nabla e\right)$ represents the changes in the macrophage-fibroblast and macrophage-macrophage adhesion strengths with respect to the growth factor and ECM in the direction of the positive gradients of growth factor and ECM, respectively. The term $- ({\bf S}_{_{mf}}f + {\bf 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 along the negative gradients of the growth factor, fibroblasts, macrophages and ECM.

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

The term $- {\bf S}_{_{me}}e\; (w_{_{g}}\nabla g + w_{_{f}}\nabla f + w_{_{m}}\nabla m + w_{_{e}}\nabla e)$ prevents the uncontrolled accumulation of macrophage-ECM, by directing them along the negative gradients of the growth factor, fibroblasts, macrophages and ECM.

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

3.   Numerical approach and simulations

Throughout this section, we investigate numerically the macrophage-ECM hypothesis mentioned above (see Cases Ⅰ and Ⅱ). 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 us denote $u: = (g, f, m, e)^{^T}$ and the right-hand side operator by $\mathrm{G}:(H^{1}(\Omega; I))^{^4}\to \mathbb{R}^{^4}$. Thus, our dynamics can now be described as {follows:

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)$ (i.e., the usual space of test functions, which coincides with $\mathcal{C}_{_{0}}^{^{\infty}}(\Omega)$ as a family of functions only) and integrate over $\Omega$. Therefore, using weak formulation, our problem is recast as follows:

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

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}({\bf x}, t) = \sum_{\tau = 0}^{l} c^{u_{i}}_{\tau}(t) \psi_{_{\tau}}\!({\bf x})$, we therefore obtain

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

where $u^{^{N}}$ is the approximation of $u(N\Delta t)$, with the uniform time discretisation step $\Delta t: = \frac{T}{N_{max}}$ and $N = 0, \dots, N_{max}$ represents 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 by taking the $L^{2}-$scalar product with respect to each basis function $\psi_{_j}$, $j = 0, \dots, l$, on both sides of the equations in system (3.3); then we apply the time discretisation outlined in system (3.5). This results in a linear system associated with the fully discretised model, as detailed in Appendix C.

3.2.   Numerical simulations

Numerical implementation. For the numerical simulations, we consider the FEM implemented in FEniCS, an open-source computing platform. In particular, here, we used 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 Subsection 3.1 and Appendix C. For details on the time-marching using the backward Euler method, see Appendix D. We ran the simulations on a time interval $\left[0, T\right]$ 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 descriptions. We assumed 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 chose the following initial conditions on the square domain $\Omega = [-1, 1]^{^{2}}$:

where ${\bf x} = (x_{_1}, x_{_2})$. Thus, the wound was assumed to be parallel to the $x_{2}$-axis; see also Figure 2.

Note that, in Appendix F, we also present 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).

3.2.1.   Numerical simulations for the local models: Gaussian kernel

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 because the cellular flux vanishes locally (see the detailed calculations in Appendix A). Figure 3 shows the spatiotemporal evolution of the growth factor $g$ (first row), fibroblasts $f$ (second row), macrophages $m$ (third row) and ECM $e$ (fourth row) at four 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 in the remodel of the ECM. Note that, at time $t = 20$, the fibroblasts and macrophages have completely invaded the wound and have 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 relative to their initial baseline levels (at $t = 0$).

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

No macrophage-ECM adhesion (Case Ⅰ: ${\bf S}_{_{me}} = 0$). First, we numerically investigate the behaviour of the local model (2.22) under the assumption that the macrophages cannot adhere to the ECM, i.e., Case Ⅰ above with ${\bf S}_{_{me}} = 0$. As before, in Figure 4, we present the spatiotemporal evolution of the growth factor $g$ (first row), fibroblasts $f$ (second row), macrophages $m$ (third row) and ECM $e$ (fourth row). In contrast to Figure 3 in which 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 that are built up in the middle of the wound ($max(f)\approx 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 throughout the 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 relative to its initial baseline level (at $t = 0$).

Including macrophage-ECM adhesion (Case Ⅱ: ${\bf S}_{_{me}} > 0$). Next, we numerically investigated the case in which 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 to remodel the ECM. Another interesting macrophage behaviour can be observed at time $t = 20$, when there are fewer macrophages inside of 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 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 ${\bf S}_{_{me}} = 0$ and Figure 5 with ${\bf 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 and macrophages: Comparison of different local models. To more clearly see the differences between various local models in terms of the fibroblasts and macrophages, in Figure 7, we plotted the densities of these two cell populations in the centre 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 the cell dynamics for (a), (b) local model (2.17); (c), (d) local model (2.22) for Case Ⅰ (${\bf S}_{_{me}} = 0$) and (e), (f) local model (2.22) for Case Ⅱ (${\bf S}_{_{me}} > 0$). In Figure 7(a), (b), we observe similar dynamics for the fibroblasts and macrophages: for $t > 5$, there are more cells in the wound region than in the surrounding tissue. This is in contrast to Figure 7(c), (d), where the fibroblast population very quickly invades 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 at around $t = 15$.

3.2.3.   Numerical simulations for the non-local models

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

Figure 8 shows the spatiotemporal evolution of the growth factor $g$ (first row), fibroblasts $f$ (second row), macrophages $m$ (third row) and ECM $e$ (fourth row) for the non-local model with the cone-shaped kernel when macrophage-ECM adhesion is allowed (Case Ⅱ). 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 to remodel the ECM. At time $t = 20$, there are many more fibroblasts inside of the wound than at the wound margins, where there are low numbers of fibroblasts.

In Figure 9, we compare Cases Ⅰ (${\bf S}_{_{me}} = 0$) and Ⅱ (${\bf S}_{_{me}} > 0$) for the non-local model (2.1) with a cone-shaped kernel (Figure 9(a)–(d)), and a Gaussian kernel (Figure 9(a')–(d')). We see that there are no significant differences between Cases Ⅰ and Ⅱ Only at $t = 100$ can we observe a slightly faster ECM remodelling process for Case Ⅰ than for Case Ⅱ (although this difference is barely noticeable, probably because ${\bf S}_{_{me}}$ is not large enough). Moreover, there are no differences between the results obtained with the 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 with the cone-shaped kernel (for Case Ⅱ: ${\bf S}_{_{me}} > 0$) at each time step, i.e., from $t = 0$ to $t = 100$. In this figure, we see that the 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 six decimal places.

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

3.2.4.   Local versus 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 Ⅱ (macrophage-ECM adhesion present) and the cone-shaped kernel.

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 that this invasion is characterised by oscillations in the ECM and fibroblast 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 passes, 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 density is still very low at the initial wound point ($x_1 = 0$), while, for the non-local model (2.1), the ECM density is much higher at $x_1 = 0$.

● 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 of the ECM densities at the wound margins (which is particularly evident at $t = 40$).

Note that in addition to the above numerical investigation of the impact of the 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 a set of simulations showing the model dynamics for $R = 0.08$, $R = 0.1$ (i.e., baseline, also investigated 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 associated with hypertrophic or keloid scars), in what follows, we numerically investigate the hypothetical case in which cells do not die due to overcrowding. In fact, in ^[55], it was observed experimentally that the fibroblasts in the centre of the keloid lesion had a reduced doubling time and lower death rates and thus achieved higher cell densities than the saturated-like densities of normal fibroblasts. Because of these observed higher fibroblast densities with reduced death rates, in what follows, 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 B):

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

● In Figure 12, we show the dynamics of the local model (2.22) (for Case Ⅱ: ${\bf S}_{_{me}} > 0$) with these new proliferation rates when we localise the cone-shaped kernel. We note that the dynamics are similar to the dynamics in Figure 5, i.e., for normal wound healing.

● In Figure 13, we show the dynamics of the non-local model (2.1) (for Case Ⅱ: ${\bf S}_{_{me}} > 0$) with the cone-shaped kernel. In contrast to Figure 12, here, we consider $e_c > 0$, which describes the possibility of cells moving down ECM gradients in the first stages of wound healing. We see that the concentration of the growth factor and the density of fibroblasts increase significantly at $t = 100$, eventually leading to the blow-up of the numerical code. The growth in the fibroblast 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 fibroblast-ECM and macrophage-ECM interactions (where only fibroblasts are assumed to contribute to ECM remodelling).

Since the non-local model (2.1) with the 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), we next investigate whether this type of fibroblast dynamics also holds for the Gaussian kernel. In Figure 14, we see that while the non-local model with the Gaussian kernel exhibits the same overgrowth of fibroblasts (for the same parameter values as in Figure 13), the corresponding local model does not show fibroblast overgrowth.

3.4.   Primary vs. secondary wound healing

As mentioned in the Introduction, wound healing can occur by primary intention (i.e., when the wound heals as the wound margins are coming together, as is the case for surgical incisions, skin grafts, or flap closures) or by secondary intention (i.e., when the wound is very large and it heals 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, 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 then at $t = 100$. To obtain these 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, the 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, ranging 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 the 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 the ECM, throughout this study we numerically investigated two possible cases: (Ⅰ) no macrophage-ECM adhesion (${\bf S}_{_{me}} = 0$) and (Ⅱ) macrophage-ECM adhesion present (${\bf S}_{_{me}} > 0$). Finally, since the wound starts to heal with the formation of the fibrin mesh, which acts as an early ECM, to describe the initial movement of cells into the wound, we chose two different adhesion functions.

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 versus non-local adhesive interactions. More precisely, the results showed that ECM remodelling is slower in the non-local models than in the local models and that ECM remodelling is slower in the presence of macrophage-ECM adhesion than when macrophage-ECM interaction 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 the production of the growth factor in specific situations: (ⅰ) when cells do not die due to overcrowding, and (ⅱ) when cells move along the 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 fibroblast 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 of these results are purely theoretical, showing the possible model behaviours in response to random variation of 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 qualitatively describe many of the observed normal and abnormal wound healing processes, particularly, 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], which can also lead to an oscillatory wave-like appearance of the tissue); the higher ECM densities sometimes observed temporarily 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 during abnormal wound healing that gives rise to hypertrophic or keloid scars and finally, the healing by primary intention (Figure 15) or 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 (which does not consider, for example, the heterogeneity of macrophages or fibroblasts inside of the normal/abnormal wounds ^[60,61]) accurately captures all aspects involved in clinical wound healing.

Future work. First, as mentioned above, we need to estimate model parameters using real 2D experimental data (e.g., using inverse problem approaches, as in ^[62]). In addition, we also need to apply data to estimate 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 biomechanical 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 results of the numerical simulations in Figure 3 (last two columns) suggest 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 the spatial 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 by using FEM. In the future it is important to stabilise these numerical schemes to ensure more accurate results (especially if we want to compare with them experimental data).

Third, 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 as 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 tools in the creation of this article.

Acknowledgments

RE and OA acknowledge funding from the French National Agency of Research (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; GR acknowledges funding from the ANR grant number ANR-21-CE45-0025-03.

Conflict of interest

The authors declare that there is no conflict of interest.

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,

Thus, for $T_{_{1}}$, we have

Following identical steps (as for $T_{_{1}}$) and changing the order of integration, we obtain that the same upper bound holds true 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:

B.   Cell growth: Logistic versus truncated logistic

Because experimental studies ^[55] have shown that keloid fibroblasts have lower death rates and reach higher densities (and, implicitly, that cell death due to overcrowding was reduced), we decided to investigate the impact of replacing the 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 the truncated logistic growth equation ($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 in which 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).

C.   FEM discretisation of non-local and local models

Using the backward (implicit) Euler method as in Eq (3.5) for the time discretisation, we have the following:

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

Next, we formulate a barycentric approximation for the non-local term $A_{_{f, m}}^{^{N}}({\bf x})$. We show only Case Ⅱ (i.e., macrophage-ECM adhesion present: $S_{_{me}} > 0$), since Case Ⅰ (i.e., no macrophage-ECM adhesion) is obtained for $S_{_{me}} = 0$.

Here, $\mathcal{A}({\bf x}): = \{T\in \mathcal{T} : \|{\bf x}-{\bf q}(T)\|_{_{2}} < R\}$ is the barycentric approximation, while ${\bf q}: \mathcal{T}\mapsto \mathbb{R}$ gives the barycentric coordinates of triangles $T\in \mathcal{T}$ i.e.,

C1.   FEM discretisation for local model (2.17), as reduced from the non-local model with the 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^{^{g}}_{_{\tau}}\}_{\tau = 1}^l, \{c^{^{f}}_{_{\tau}}\}_{\tau = 1}^l, \{c^{^{m}}_{_{\tau}}\}_{\tau = 1}^l, \{c^{^{e}}_{_{\tau}}\}_{\tau = 1}^l$ such that

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

represents their gradients, where

$S^{^{\kappa}}_{_{g}} = \frac{\partial}{\partial g} S_{\kappa} = S^{^{\kappa}}_{_{e}} = \frac{\partial}{\partial e} S_{\kappa} = \frac{1}{\varepsilon^2}$ and $\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)$ discretised in time and space, i.e.,

C2.   FEM discretisation for local model (2.22), as reduced from non-local model with the cone-shaped kernel

Below, we show the FEM discretisation for the local models. However, since Case Ⅰ (i.e., no macrophage-ECM adhesion: $S_{_{me}} = 0$) is just a simplification of Case Ⅱ (which includes macrophage-ECM adhesion), in what follows, we present only the FEM discretisation for Case Ⅱ ($S_{_{me}} > 0$).

The weak form of system (2.22) as 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^{^{g}}_{_{\tau}}\}_{\tau = 1}^l, \{c^{^{f}}_{_{\tau}}\}_{\tau = 1}^l, \{c^{^{m}}_{_{\tau}}\}_{\tau = 1}^l, \{c^{^{e}}_{_{\tau}}\}_{\tau = 1}^l$ such that

D.   Time-marching in FEniCS

We rewrite Eq (3.5) in the standard form $F\left(u^{^{N\!+\!1}}; u^{^{N}}; v\right) = 0,$ where

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

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

Set some stopping time $T$

$t = \Delta t$

while $t\leq 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\xleftarrow{} t+ \Delta t$ (update time)

  $z_{_{prev}}\xleftarrow{}z$ (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 by using the fine mesh $256\times256$ for our exact solution while mesh sizes $4\times 4, \; 8\times 8, \; 16\times 16, \; 32\times 32$, $64\times 64,$ and $128\times 128$ represented our approximate solutions. Figure A2 is a plot of the L$_2$ errors obtained against the different mesh sizes. The errors observed using mesh sizes $32\times 32$ are insignificant compared to those of mesh sizes $4\times 4$ and $8\times 8$. For this reason, we decided to use a mesh size of $32\times 32$ for all simulations presented in this study.

F.   Simulation results for two new initial conditions

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

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

In Figures A3 and A4, we see that at time $t = 100$, both of the wounds heal normally as the ECM density approaches its maximum level while the fibroblast and macrophage tissue levels approach their baseline levels. However, unlike the healing in the previous figures (i.e., for the original initial conditions), here, we see that, just before the final healing, the ECM density is increased in the wound region (i.e., above the density of the ECM density in the surrounding tissue), and then it decreases towards its density in the surrounding tissue. We note that such transient behaviour could have been possible for the simulations performed under the original initial conditions. However, since we did not show the solutions at every single time step, such detailed behaviour was lost before and was randomly revealed through Figures A3 and A4. This suggests that a more thorough analytical and numerical investigation of the transient model dynamics is necessary to understand the behaviour of these non-local mathematical models.

G.   Sensitivity test for the sensing radius

Since there are some small differences between the non-local models and the corresponding local models (i.e., as obtained in the limit $R\to 0$), in what follows, we present a numerical sensitivity test for the shape of the solutions for different values of the cell sensing radius $R$. In Figure A5, we show a series of space slices 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 remodelling of the ECM and a faster decay of the fibroblasts, macrophages and the growth factor.