# Upscaling between an Agent-Based Model (Smoothed Particle Approach) and a
Continuum-Based Model for Skin Contractions

Q. Peng, F.J. Vermolen

###### Abstract

Skin contraction is an important biophysical process that takes place during
and after recovery of deep tissue injury. This process is mainly caused by
fibroblasts (skin cells) and myofibroblasts (differentiated fibroblasts) that
exert pulling forces on the surrounding extracellular matrix (ECM). Modelling
is done in multiple scales: agent–based modelling on the microscale and
continuum–based modelling on the macroscale. In this manuscript we present
some results from our study of the connection between these scales. For the
one–dimensional case, we managed to rigorously establish the link between the
two modelling approaches for both closed–form solutions and finite–element
approximations. For the multi–dimensional case, we computationally evidence
the connection between the agent–based and continuum–based modelling
approaches.

Key Words: traction forces, mechanics, smoothed particle approach, agent-based
modelling, finite element methods, continuum-based modelling

##  1 Introduction

Wound healing is a spontaneous process for the skin to cure itself after an
injury. It is a complicated combination of various cellular events that
contribute to resurfacing, reconstituting and restoring of the tensile
strength of the injured skin. For superficial wounds that only stay in
epidermis, the wound can be healed without any problem. However, for severe
injuries, in particular dermal wound, they may result in various pathological
problems, such as contractures.

Contractures concur with disabilities and disfunctioning of patients, which
cause a significantly severe impact on patients’ daily life. Contractures are
recognized as excessive and problematic contractions, which occur due to the
pulling forces exerted by the (myo)fibroblasts on their direct environment,
i.e. the extra cellular matrix (ECM). Contractions mainly start occurring in
the proliferation phase of wound healing: the fibroblasts are migrating
towards the wound and differentiating into myofibroblasts due to the high
concentration of transforming growth factor-beta (TGF-beta). Usually, $5-10\%$
reduction of wound area has been observed in clinical trials. A more detailed
biological description can be found in [5, 4, 7, 9].

In our previous work [13], a formalism to describe the mechanism of the
displacement of the ECM has been used, which is firstly developed by Boon et
al. [2] and improved further by Koppenol [8]. Regarding the elasticity
equation with point forces, we realized that the solution to the partial
differential equation is singular for dimensionality exceeding one. Hence, we
developed various alternatives to improve the accuracy of the solution in [11,
12].

We have been working with agent-based models so far, which model the cells as
individuals and define the formalism of pulling forces by superposition
theory. However, once the wound scale is larger, the agent-based model is
increasingly expensive from a computational perspective, and hence, the cell
density model is preferred, which considers many cells as one collection in a
unit. In this manuscript, we investigate and discover the connections between
these two models, in the perspective of modelling the mechanism of pulling
forces exerted by the (myo)fibroblasts. As the consistency between the
smoothed particle approach (SPH approach) and the immersed boundary approach
has been proven both analytically and numerically [11, 12], we select the SPH
approach here due to its continuity and smoothness, to compare with the cell
density model using finite-element methods.

The manuscript is structured as follows. We start introducing both models in
one dimension in Section 2, then in Section 3 we extend the models to two
dimensions. Section 4 displays the numerical results in one and two
dimensions. Finally, some conclusions are shown in Section 5.

##  2 Mathematical Models in One Dimension

Considering one-dimensional force equilibrium, the equations are given by

$\displaystyle-\frac{d\sigma}{dx}$ $\displaystyle=f,$ Equation of Equlibirum,
$\displaystyle\epsilon$ $\displaystyle=\frac{du}{dx},$ Strain-Displacement
Relation, $\displaystyle\sigma$ $\displaystyle=E\epsilon,$ Constitutive
Equation.

By substituting $E=1$, the equations above can be combined to Laplacian
equation in one dimension:

$-\frac{d^{2}u}{dx^{2}}=f.$ (2.1)

###  2.1 Smoothed Particle Approach

In [12], a smoothed particle approach (SPH approach) is developed as an
alternative of the Dirac Delta distribution describing the point forces
exerted by the biological cells, in the application of wound healing:

$(BVP_{SPH})\left\\{\begin{aligned}
-\frac{d^{2}u}{dx^{2}}&=P_{SPH}\sum_{i=1}^{N_{s}}\delta^{\prime}_{\varepsilon}(x-s_{i}),\mbox{$x\in(0,L)$,}\\\
u(0)&=u(L)=0,\end{aligned}\right.$ (2.2)

where $P_{SPH}$ is the magnitude of the forces, $\delta_{\varepsilon}(x)$ is
the Gaussian distribution with variance $\varepsilon$ and $s_{i}$ is the
centre position of biological cell $i$. One can solve the partial differential
equations (PDEs) with finite-element methods. The corresponding weak form is
given by

$(WF\textsubscript{SPH})\left\\{\begin{aligned} &\text{Find $u\in
H_{0}^{1}((0,L))$, such that}\\\
&\int_{0}^{L}u^{\prime}\phi^{\prime}dx=\int_{0}^{L}\sum_{i=1}^{N_{s}}P_{SPH}\delta^{\prime}_{\varepsilon}(x-s_{i})\phi
dx,\\\ &\text{for all $\phi\in H_{0}^{1}((0,L))$.}\end{aligned}\right.$

Without this knowledge, the existence and uniqueness of the
$H^{1}_{0}$-solution follows as well from the application of the Lax–Milgram
theorem [3], where it is immediately obvious that the bilinear form in the
left–hand side is symmetric and positive definite.

###  2.2 Cell Density Approach

A cell density approach is often used in the large scale, so that the
computational efficiency is much improved compared with the agent-based model.
According to the model in [8], the force in two dimensions can be determined
by the divergence of $n_{c}\cdot\bm{I}$, where $n_{c}$ is the local density of
the biological cells and $\bm{I}$ is the identity tensor. In one dimension,
the cell density approach is expressed as:

$(BVP_{den})\left\\{\begin{aligned}
-\frac{d^{2}u}{dx^{2}}&=P_{den}\frac{dn_{c}}{dx},\mbox{$x\in(0,L)$,}\\\
u(0)&=u(L)=0,\end{aligned}\right.$ (2.3)

where $P_{den}$ is the magnitude of the forces. The corresponding weak form is
given by

$(WF\textsubscript{den})\left\\{\begin{aligned} &\text{Find $u\in
H_{0}^{1}((0,L))$, such that}\\\
&\int_{0}^{L}u^{\prime}\phi^{\prime}dx=\int_{0}^{L}P_{den}n^{\prime}_{s}\phi
dx,\\\ &\text{for all $\phi\in H_{0}^{1}((0,L))$.}\end{aligned}\right.$

###  2.3 Consistency between Two Models

####  2.3.1 Analytical Solutions with Specific Locations of Biological Cells

To express the analytical solution, it is necessary to determine the locations
of the biological cells, such that the cell density can be written as an
analytical function of the positions. We assume, there are $N_{s}$ cells
distributed uniformly in the subdomain $(a,b)$ of the computational domain
$(0,L)$. Hence, the distance between the center position of any two adjacent
biological cells is constant, which we denote $\Delta s=(b-a)/N_{s}$ and the
first and the $N_{s}$-th cell are located at $x=a+\Delta s/2$ and $x=b-\Delta
s/2$, respectively. With homogeneous Dirichlet boundary conditions, and
suppose $P_{SPH}=P\Delta s$ and variance $\varepsilon=\Delta s$, the boundary
value problem of the SPH approach is expressed as

$(BVP^{1}_{SPH})\left\\{\begin{aligned} -\frac{d^{2}u_{1}}{dx^{2}}&=P\Delta
s\sum_{i=1}^{N_{s}}\delta^{\prime}_{\Delta s}(x-s_{i}),\mbox{$x\in(0,L)$,}\\\
u_{1}(0)&=u_{1}(L)=0,\end{aligned}\right.$ (2.4)

where $P$ is a positive constant and $s_{i}$ is the centre position of the
biological cells. Utilizing the superposition principle, the analytical
solution is given by

$u_{1}(x)=P\Delta
s\sum_{i=1}^{N_{s}}\frac{1}{2}\left\\{\left(\frac{x}{L}-1\right)\operatorname{erf}\left(\frac{s_{i}}{\sqrt{2}\Delta
s}\right)+\frac{x}{L}\operatorname{erf}\left(\frac{L-s_{i}}{\sqrt{2}\Delta
s}\right)-\operatorname{erf}\left(\frac{x-s_{i}}{\sqrt{2}\Delta
s}\right)\right\\},$ (2.5)

where $\operatorname{erf}(x)$ is the error function defined as
$\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp(-t^{2})dt$ [15].
Since the biological cells are uniformly located between $a$ and $b$
($0<a<b<L$), $\frac{dn_{c}}{dx}$ can be rephrased as

$\frac{dn_{c}}{dx}=\left\\{\begin{aligned}
&\frac{1}{t},&a-\frac{t}{2}<x<a+\frac{t}{2},\\\
&-\frac{1}{t},&b-\frac{t}{2}<x<b+\frac{t}{2},\\\
&0,&\mbox{otherwise,}\end{aligned}\right.$

where $t$ is a small positive constant. Taking $t$ to zero, the above
expression converges to $\delta(x-a)-\delta(x-b)$. Hence, the boundary value
problem of the cell density model can be written as

$(BVP^{1}_{den})\left\\{\begin{aligned}
-\frac{d^{2}u_{2}}{dx^{2}}&=P\frac{dn_{c}}{dx}\rightarrow
P(\delta(x-a)-\delta(x-b)),\mbox{$x\in(0,L)$,}\\\
u_{2}(0)&=u_{2}(L)=0,\end{aligned}\right.$ (2.6)

where $\delta(x)$ is the Dirac Delta distribution and $a$ and $b$ are the left
and right endpoint of the subdomain (where biological cells are uniformly
located) respectively. The analytical solution is then expressed as

$u_{2}(x)=P(G(x,a)-G(x,b)),$ (2.7)

where $G(x,x^{\prime})$ is the Green’s function [6], defined by

$G(x,x^{\prime})=(1-\frac{x^{\prime}}{L})x-max(x-x^{\prime},0),$

in the computational domain $(0,L)$.

Actually, the convergence between $u_{1}(x)$ and $u_{2}(x)$ can be proven as
$\Delta s\rightarrow 0^{+}$ by a proposition. Firstly, we introduce
Chebyshev’s Inequality:

######  Lemma 2.1.

(Chebyshev’s Inequality [10]) Denote $X$ as a random variable with finite mean
$\mu$ and finite variance $\sigma^{2}$. Then for any positive
$k\in\mathbb{R}$, the following inequality holds:

$\mathbb{P}(|X-\mu|\geqslant k)\leqslant\frac{\sigma^{2}}{k^{2}},$

where $\mathbb{P}(A)$ is the probability of event $A$. The above inequality
can also be rephrased as

$\mathbb{P}(|X-\mu|\leqslant k)\geqslant 1-\frac{\sigma^{2}}{k^{2}}.$

######  Proposition 2.1.

Let $u_{1}(x)$ as described in Eq (2.4) be the exact solution to
$(BVP^{1}_{SPH})$ and $u_{2}(x)$ as described in Eq (2.6) be the exact
solution to $(BVP_{den}^{1})$. As $\Delta s\rightarrow 0^{+}$, $u_{1}(x)$
converges to $u_{2}(x)$.

###### Proof.

For the standard Gaussian distribution in one dimension, the cumulative
distribution function is given by $\displaystyle
F(x)=\frac{1}{2}\left(1+\operatorname{erf}\left(\frac{x}{\sqrt{2}}\right)\right)$.
Thus, we obtain

$\operatorname{erf}\left(\frac{x}{\sqrt{2}}\right)=2F(x)-1.$ (2.8)

By Chebyshev’s Inequality (see Lemma 2.1), one can conclude that for any
positive $k$,

$F(k)-F(-k)\geqslant 1-\frac{1}{k^{2}}.$ (2.9)

Note that $1-F(k)=F(-k)$ due to the symmetry of standard Gaussian
distribution. Hence, Eq (2.9) implies

$\displaystyle F(k)\geqslant 1-\frac{1}{k^{2}}+F(-k)=1-\frac{1}{k^{2}}+1-F(k)$
$\displaystyle\Leftrightarrow$ $\displaystyle 1-\frac{1}{k^{2}}\leqslant
F(k)\leqslant 1,$

and analogously, $0\leqslant F(-k)\leqslant\frac{1}{k^{2}}$ is implied.
Together with Eq (2.8), it gives

$\left\\{\begin{aligned}
&1-\frac{1}{k^{2}}\leqslant\operatorname{erf}\left(\frac{k}{\sqrt{2}}\right)\leqslant
1,\\\
&-1\leqslant\operatorname{erf}\left(-\frac{k}{\sqrt{2}}\right)\leqslant-1+\frac{1}{k^{2}}.\end{aligned}\right.$

Let $k=\frac{s_{i}}{\Delta s}>0$, for any
$s_{i}\in(a,b)\subset(0,L),i=\\{1,\dots,N_{s}\\}$, then

$\left\\{\begin{aligned} &1-\left(\frac{\Delta
s}{s_{i}}\right)^{2}\leqslant\operatorname{erf}\left(\frac{s_{i}}{\sqrt{2}\Delta
s}\right)\leqslant 1,\\\
&-1\leqslant\operatorname{erf}\left(-\frac{s_{i}}{\sqrt{2}\Delta
s}\right)\leqslant-1+\left(\frac{\Delta
s}{s_{i}}\right)^{2}.\end{aligned}\right.$

As it has been defined earlier that $\Delta s=(b-a)/N_{s}$, we obtain

$\displaystyle\sum_{i=1}^{N_{s}}\left(1-\left(\frac{\Delta
s}{s_{i}}\right)^{2}\right)\Delta
s\leqslant\sum_{i=1}^{N_{s}}\operatorname{erf}\left(\frac{s_{i}}{\sqrt{2}\Delta
s}\right)\Delta s\leqslant\sum_{i=1}^{N_{s}}\Delta s$
$\displaystyle\Rightarrow$ $\displaystyle(b-a)-(\Delta
s)^{3}\sum_{i=1}^{N_{s}}\frac{1}{s_{i}^{2}}\leqslant\sum_{i=1}^{N_{s}}\operatorname{erf}\left(\frac{s_{i}}{\sqrt{2}\Delta
s}\right)\Delta s\leqslant(b-a).$

Since $\displaystyle\lim_{\Delta s\rightarrow 0^{+}}(\Delta
s)^{3}\sum_{i=1}^{N_{s}}\frac{1}{s_{i}^{2}}=0$ for any
$s_{i}\in(a,b)\subset(0,L)$, the Squeeze Theorem [1] implies that

$\displaystyle\lim_{\Delta s\rightarrow
0^{+}}\sum_{i=1}^{N_{s}}\operatorname{erf}\left(\frac{s_{i}}{\sqrt{2}\Delta
s}\right)\Delta s=b-a.$ (2.10)

Analogously, we obtain that for any $s_{i}\in(a,b)\subset(0,L)$,

$\displaystyle\lim_{\Delta s\rightarrow
0^{+}}\sum_{i=1}^{N_{s}}\operatorname{erf}\left(-\frac{s_{i}}{\sqrt{2}\Delta
s}\right)\Delta s=a-b.$ (2.11)

Thus, it can be concluded that for any series of real number
$\\{x_{i}\\}\in\mathbb{R}^{n}$, when $x_{i+1}-x_{i}=\Delta s$ and $x_{i}$ is
either all positive or all negative for any $i=\\{1,\cdots,N_{s}\\}$,

$\displaystyle\lim_{\Delta s\rightarrow
0^{+}}\sum_{i=1}^{N_{s}}\operatorname{erf}\left(\frac{x_{i}}{\sqrt{2}\Delta
s}\right)\Delta s=(b-a)\operatorname{sgn}(x_{i}),$ (2.12)

where $\operatorname{sgn}(x)$ is sign function defined by

$\operatorname{sgn}(x)=\left\\{\begin{aligned} &1,\mbox{ if $x>0$,}\\\
&0,\mbox{ if $x=0$,}\\\ &-1,\mbox{ if $x<0$.}\\\ \end{aligned}\right.$

We rewrite $u_{1}(x)$ as

$\displaystyle
u_{1}(x)=P\left[\frac{1}{2}\left(\frac{x}{L}-1\right)\sum_{i=1}^{N_{s}}\operatorname{erf}\left(\frac{s_{i}}{\sqrt{2}\Delta
s}\right)\Delta
s+\frac{1}{2}\frac{x}{L}\sum_{i=1}^{N_{s}}\operatorname{erf}\left(\frac{L-s_{i}}{\sqrt{2}\Delta
s}\right)\Delta
s-\frac{1}{2}\sum_{i=1}^{N_{s}}\operatorname{erf}\left(\frac{x-s_{i}}{\sqrt{2}\Delta
s}\right)\Delta s\right].$

Combining Eq (2.10), (2.11) and (2.12), $u_{1}(x)$ is given by

$\displaystyle u_{1}(x)$
$\displaystyle=P\left[\frac{1}{2}\left(\frac{x}{L}-1\right)(b-a)+\frac{1}{2}\frac{x}{L}(b-a)+\frac{1}{2}[(x-b)\operatorname{sgn}(x-b)-(x-a)\operatorname{sgn}(x-a)]\right]$
$\displaystyle=P\left[\left(\frac{x}{L}-\frac{1}{2}\right)(b-a)+\frac{1}{2}[(x-b)\operatorname{sgn}(x-b)-(x-a)\operatorname{sgn}(x-a)]\right]$
$\displaystyle=P\left[\left(\frac{x}{L}-\frac{1}{2}\right)(b-a)+\frac{1}{2}[|x-b|-|x-a|]\right]$
$\displaystyle=\left\\{\begin{aligned} &P\frac{x}{L}(b-a),\quad 0\leqslant
x\leqslant a,\\\ &P\frac{x}{L}(b-a)-x+a,\quad a<x\leqslant b,\\\
&P\left(\frac{x}{L}-1\right)(b-a),\quad b<x\leqslant L,\\\
\end{aligned}\right.$

Rewriting $u_{2}(x)$ regarding different domain gives

$u_{2}(x)=\left\\{\begin{aligned} &P\frac{x}{L}(b-a),\quad 0\leqslant
x\leqslant a,\\\ &P\frac{x}{L}(b-a)-x+a,\quad a<x\leqslant b,\\\
&P\left(\frac{x}{L}-1\right)(b-a),\quad b<x\leqslant L.\\\
\end{aligned}\right.$

Hence, we conclude that $u_{1}(x)$ converges to $u_{2}(x)$ as $\Delta
s\rightarrow 0^{+}$. ∎

####  2.3.2 Finite-Element Method Solutions with Arbitrary Locations of
Biological Cells

For the finite-element method, we select the piecewise Lagrangian linear basis
functions. We divide the computational domain into $N_{e}$ mesh elements, with
the nodal point $x_{1}=0$ and $x_{N_{e}+1}=L$. For the implementation, we
define the cell density as the count of biological cell in every mesh element
divided by the length of the mesh element, hence, it is a constant within
every mesh element. In other words, in the mesh element $[x_{j},x_{j+1}]$, the
count of the biological cell is defined by

$N_{c}([x_{j},x_{j+1}])=\int_{x_{j}}^{x_{j+1}}n_{c}([x_{j},x_{j+1}])dx=hn_{c}([x_{j},x_{j+1}]),$

for any $j\in\\{1,\dots,N_{e}\\}$, where $h$ is the size of every mesh
element. Different from $(BVP^{1}_{SPH})$ where $\Delta s$ is the variance of
$\delta_{\varepsilon}$, for finite-element methods, we set $\varepsilon=h/3$,
such that the integration of $\delta_{h/3}(x-x^{\prime})$ for any
$0<x^{\prime}<L$ over any mesh element with size $h$, is close to $1$ (see
Lemma 2.2). With the two approaches, the boundary value problems with
Dirichlet boundary condition are defined by

$(BVP^{2}_{SPH})\left\\{\begin{aligned}
-\frac{d^{2}u_{1}}{dx^{2}}&=Ph\sum_{i=1}^{N_{s}}\delta^{\prime}_{h/3}(x-s_{i}),\mbox{$x\in(0,L)$,}\\\
u_{1}(0)&=u_{1}(L)=0,\end{aligned}\right.$ (2.13)

and

$(BVP^{2}_{den})\left\\{\begin{aligned}
-\frac{d^{2}u_{2}}{dx^{2}}&=Ph\frac{dn_{c}}{dx},\mbox{$x\in(0,L)$,}\\\
u_{2}(0)&=u_{2}(L)=0,\end{aligned}\right.$ (2.14)

where $s_{i}$ is the position of biological cells, $h$ is the mesh size and
$N_{s}$ is the total number of cells in the computational domain. The
consistency between $(BVP^{2}_{SPH})$ and $(BVP^{2}_{den})$ can be verified by
the following lemma and theorem.

######  Lemma 2.2.

(Empirical rule [14]) Given the Gaussian distribution of mean $\mu$ and
variance $\varepsilon$:

$\delta_{\varepsilon}(x-\mu)=1/\sqrt{2\pi\varepsilon^{2}}\exp\\{-(x-\mu)^{2}/(2\varepsilon^{2})\\},$

then the following integration can be computed:

  1. 1.

$\int_{\mu-\varepsilon}^{\mu+\varepsilon}\delta_{\varepsilon}(x-\mu)dx\approx
0.6827;$

  2. 2.

$\int_{\mu-2\varepsilon}^{\mu+2\varepsilon}\delta_{\varepsilon}(x-\mu)dx\approx
0.9545;$

  3. 3.

$\int_{\mu-3\varepsilon}^{\mu+3\varepsilon}\delta_{\varepsilon}(x-\mu)dx\approx
0.9973.$

######  Theorem 2.1.

Denote $u^{h}_{1}(x)$ and $u^{h}_{2}(x)$ respectively the solution to
$(BVP^{2}_{SPH})$ and $(BVP^{2}_{den})$. With Lagrangian linear basis
functions for the finite element method, $u^{h}_{1}(x)$ converges to
$u^{h}_{2}(x)$, as the size of the mesh element $h\rightarrow 0^{+}$,
regardless of the positions of biological cells.

###### Proof.

We define $v^{h}(x)=u^{h}_{1}(x)-u^{h}_{2}(x)$, then the boundary value
problem to solve $v^{h}(x)$ is given by

$(BVP^{1}_{v})\left\\{\begin{aligned}
-\frac{d^{2}v^{h}}{dx^{2}}&=Ph\left(\sum_{i=1}^{N_{s}}\delta^{\prime}_{h/3}(x-s_{i})-\frac{dn_{c}}{dx}\right),\mbox{$x\in(0,L)$,}\\\
v(0)&=v(L)=0.\end{aligned}\right.$ (2.15)

The corresponding Galerkin’s form reads as

$(GF_{v}^{1})\left\\{\begin{aligned} &\text{Find $v^{h}\in H_{0}^{1}((0,L))$,
such that}\\\
&\int_{0}^{L}\frac{dv^{h}}{dx}\phi^{\prime}dx=\int_{0}^{L}Ph\left(\sum_{i=1}^{N_{s}}\delta^{\prime}_{h/3}(x-s_{i})-\frac{dn_{c}}{dx}\right)\phi
dx,\\\ &\text{for all $\phi\in H_{0}^{1}((0,L))$.}\end{aligned}\right.$

Using integration by parts and letting
$\phi=\phi_{j},j\in\\{1,\dots,N_{e}+1\\}$, the equation in $(GF_{v}^{1})$ can
be rewritten by

$\displaystyle\int_{0}^{L}\frac{dv^{h}}{dx}\phi^{\prime}_{j}dx$
$\displaystyle=\int_{0}^{L}Ph\left(\sum_{i=1}^{N_{s}}\delta^{\prime}_{h/3}(x-s_{i})-\frac{dn_{c}}{dx}\right)\phi_{j}dx$
$\displaystyle=\left[Ph\sum_{i=1}^{N_{s}}\delta_{h/3}(x-s_{i})\phi_{j}\right]^{L}_{0}-[Phn_{c}\phi_{j}]^{L}_{0}-\int_{0}^{L}Ph\left(\sum_{i=1}^{N_{s}}\delta_{h/3}(x-s_{i})-n_{c}\right)\phi^{\prime}_{j}dx$
(Boundary condition)
$\displaystyle=-\int_{0}^{L}Ph\left(\sum_{i=1}^{N_{s}}\delta_{h/3}(x-s_{i})-n_{c}\right)\phi^{\prime}_{j}dx$
$\displaystyle=Ph\sum_{j=1}^{N_{e}}\left\\{\int_{x_{j}}^{x_{j+1}}n_{c}\phi_{j}^{\prime}dx-\int_{x_{j}}^{x_{j+1}}\sum_{i=1}^{N_{s}}\delta_{h/3}(x-s_{i})\phi^{\prime}_{j}dx\right\\}$
$\displaystyle=P\sum_{j=1}^{N_{e}}\left[N_{c}([x_{j},x_{j+1}])-\int_{x_{j}}^{x_{j+1}}\sum_{i=1}^{N_{s}}\delta_{h/3}(x-s_{i})dx\right]$
($h\rightarrow 0^{+}$, Lemma 2.2) $\displaystyle\rightarrow 0,$

since it can be defined that $\displaystyle
N([x_{j},x_{j+1}])=\int_{x_{j}}^{x_{j+1}}\sum_{i=1}^{N_{s}}\delta(x-s_{i})dx$.
∎

##  3 Mathematical Models in Two Dimensions

###  3.1 Smoothed Particle Approach and Cell Density Approach

In multi dimensional case, the equation of conservation of momentum over the
computational domain $\Omega$, without considering inertia, is given by

$-\nabla\cdot\bm{\sigma}=\bm{f}.$

We consider a linear, homogeneous and isotropic domain, with Hooke’s Law, the
stress tensor $\bm{\sigma}$ is defined as

$\bm{\sigma}=\frac{E}{1+\nu}\left\\{\bm{\epsilon}+\operatorname{tr}(\epsilon)\left[\frac{\nu}{1-2\nu}\right]\bm{I}\right\\},$
(3.1)

where $E$ is the Young’s modulus of the material, $\nu$ is Poisson’s ratio and
$\bm{\epsilon}$ is the infinitesimal strain tensor:

$\bm{\epsilon}=\frac{1}{2}[\bm{\nabla u}+(\bm{\nabla u})^{T}].$

Considering a subdomain $\Omega_{w}\subset\Omega$, where the center positions
of the biological cells are located, then the SPH approach and cell density
approach with homogeneous Dirichlet boundary condition are derived by

$(BVP^{3}_{SPH})\left\\{\begin{aligned}
-\nabla\cdot\bm{\sigma}&=P_{SPH}\sum_{i=1}^{N_{s}}\nabla\delta_{\varepsilon}(\bm{x}-\bm{s_{i}}),\mbox{$\bm{x}\in\Omega$,}\\\
\bm{u_{1}}(\bm{x})&=\bm{0},\mbox{$\bm{x}\in\partial\Omega$,}\end{aligned}\right.$
(3.2)

and

$(BVP^{3}_{den})\left\\{\begin{aligned}
-\nabla\cdot\bm{\sigma}&=P_{den}\nabla\cdot(n_{c}\bm{I}),\mbox{$\bm{x}\in\Omega$,}\\\
\bm{u_{2}}(\bm{x})&=\bm{0},\mbox{$\bm{x}\in\partial\Omega$.}\end{aligned}\right.$
(3.3)

###  3.2 Consistency between Two Approaches in Finite-Element Method

To prove the consistency between these two approaches, we define that for the
triangular mesh element $e_{k},k\in\\{1,\dots,N_{e}\\}$, where $N_{e}$ is the
total number of mesh elements in $\Omega$, the density of biological cell
$n_{c}(e_{k})$ is constant and the count of biological cells $N_{c}(e_{k})$ is
expressed by

$N_{c}(e_{k})=\int_{e_{k}}n_{c}(e_{k})d\Omega=A(e_{k})n_{c}(e_{k}),$

where $A(e_{k})$ is the area of mesh element $e_{k}$. Similar to the process
of proof in one dimension, we state the following theorem:

######  Theorem 3.1.

Denote $\bm{u^{h}_{1}}(\bm{x})$ and $\bm{u^{h}_{2}}(\bm{x})$ respectively the
solution to $(BVP^{3}_{SPH})$ and $(BVP^{3}_{den})$ with $P_{SPH}=P_{den}=P$.
With Lagrange linear basis functions for the finite element method,
$\bm{u^{h}_{1}}(\bm{x})$ converges to $\bm{u^{h}_{2}}(\bm{x})$, as the size of
a triangular mesh element is taken to zero, regardless of the positions of
biological cells.

###### Proof.

We consider
$\bm{v^{h}}(\bm{x})=\bm{u^{h}_{1}}(\bm{x})-\bm{u^{h}_{2}}(\bm{x})$, then the
boundary value problem to solve $\bm{v^{h}}(\bm{x})$ is given by

$(BVP^{3}_{v})\left\\{\begin{aligned}
-\nabla\cdot\bm{\sigma}&=P\left[\sum_{i=1}^{N_{s}}\nabla\delta_{\varepsilon}(\bm{x}-\bm{s_{i}})-\nabla\cdot(n_{c}\bm{I})\right],\mbox{$\bm{x}\in\Omega$,}\\\
\bm{v}(\bm{x})&=\bm{0},\mbox{$\bm{x}\in\partial\Omega$}.\end{aligned}\right.$
(3.4)

The corresponding Galerkin’s form reads as

$(GF^{3}_{v})\left\\{\begin{aligned} &\text{Find
$\bm{v^{h}}\in\bm{H_{0}^{1}}(\Omega)$, such that}\\\
&\int_{\Omega}\bm{\sigma}(\bm{v^{h}}):\nabla\bm{\phi^{h}}d\Omega=\int_{\Omega}P\left[\sum_{i=1}^{N_{s}}\nabla\delta_{h}(\bm{x}-\bm{s_{i}})-\nabla\cdot(n_{c}\bm{I})\right]\cdot\bm{\phi^{h}}d\Omega,\\\
&\text{for all $\bm{\phi^{h}}\in\bm{H_{0}^{1}}(\Omega)$.}\end{aligned}\right.$

With integral by parts and letting
$\bm{\phi^{h}}=\bm{\phi^{h}_{k}},k\in\\{1,\dots,N_{e}\\}$, the equation in
$(GF_{v}^{3})$ can be rewritten by

$\displaystyle\int_{\Omega}\bm{\sigma}(\bm{v^{h}})\nabla\bm{\phi^{h}_{k}}d\Omega$
$\displaystyle=\int_{\Omega}P\left[\sum_{i=1}^{N_{s}}\nabla\delta_{\varepsilon}(\bm{x}-\bm{s_{i}})-\nabla\cdot(n_{c}\bm{I})\right]\cdot\bm{\phi^{h}_{k}}d\Omega$
$\displaystyle=P\left\\{\left[\int_{\partial\Omega}\sum_{i=1}^{N_{s}}\delta_{\varepsilon}(\bm{x}-\bm{s_{i}})\bm{\phi^{h}_{k}}\bm{n}d\Gamma-\int_{\Omega}\sum_{i=1}^{N_{s}}\delta_{\varepsilon}(\bm{x}-\bm{s_{i}})\nabla\cdot\bm{\phi^{h}_{k}}d\Omega\right]\right.$
$\displaystyle\left.-\left[\int_{\partial\Omega}n_{c}(\bm{I\phi^{h}_{k}})\cdot\bm{n}d\Gamma-\int_{\Omega}n_{c}\bm{I}:\nabla\bm{\phi^{h}_{k}}d\Omega\right]\right\\}$
(Boundary condition)
$\displaystyle=-P\int_{\Omega}\sum_{i=1}^{N_{s}}\delta_{\varepsilon}(\bm{x}-\bm{s_{i}})\nabla\cdot\bm{\phi^{h}_{k}}-n_{c}\bm{I}:\nabla\bm{\phi^{h}_{k}}d\Omega$
$\displaystyle=-P\sum_{k=1}^{N_{e}}\int_{e_{k}}\sum_{i=1}^{N_{s}}\delta_{\varepsilon}(\bm{x}-\bm{s_{i}})\nabla\cdot\bm{\phi^{h}_{k}}-n_{c}\bm{I}:\nabla\bm{\phi^{h}_{k}}d\Omega$
$\displaystyle=-P\sum_{k=1}^{N_{e}}\nabla\cdot\bm{\phi^{h}_{k}}\int_{e_{k}}\sum_{i=1}^{N_{s}}\delta_{\varepsilon}(\bm{x}-\bm{s_{i}})-n_{c}d\Omega$
$\displaystyle=P\sum_{k=1}^{N_{e}}\nabla\cdot\bm{\phi^{h}_{k}}\left[N_{c}(e_{k})-\int_{e_{k}}\sum_{i=1}^{N_{s}}\delta_{\varepsilon}(\bm{x}-\bm{s_{i}})d\Omega\right]$
($\varepsilon\rightarrow 0^{+}$) $\displaystyle\rightarrow 0,$

since it can be defined that $\displaystyle
N_{c}(e_{k})=\int_{e_{k}}\sum_{i=1}^{N_{s}}\delta(\bm{x}-\bm{s_{i}})d\Omega$.
Note that in two dimensions, $\varepsilon$ needs to be sufficiently small
compared to the size of a triangular mesh element.

∎

##  4 Results

Results in both one and two dimensions are discussed in this section. Since
the objective of this manuscript is to investigate the consistency and the
connections between the SPH approach and the cell density approach, all the
parameters are dimensionless.

###  4.1 One-Dimensional Results

We show the results by analytical solutions in Figure 4.1 with various values
of $\Delta s$ (i.e. depending on different number of biological cells in the
subdomain $(a,b)$). Here, the computational domain is $(0,7)$ with $L=7$ and
the subdomain where the biological cells locate uniformly is $(2,5)$ with
$a=2$ and $b=5$. With the decrease of the variance in the Gaussian
distribution in $(BVP^{2}_{SPH})$, the curves gradually overlap, which
verifies the convergence between the analytical solutions to these two
approaches.

(a) $\Delta s=0.3$

(b) $\Delta s=0.06$

(c) $\Delta s=0.03$

(d) $\Delta s=0.006$

Figure 4.1: The exact solutions to $(BVP^{1}_{SPH})$ and $(BVP^{1}_{den})$ are
shown, with various values of $\Delta s$, which is the distance between centre
positions of any two adjacent biological cells. Blue points are the centre
positions of biological cells. Red curves represent the solutions to
$(BVP^{1}_{SPH})$ and blue curves represent the solutions to
$(BVP^{1}_{den})$.

To implement the model, there are two different algorithms shown in Figure 4.2
and 4.3. Depending on different circumstances, the implementation method is
elected. The cell density in one dimension is defined as the number of cells
per length unit. In other words, the cell count in a given domain can be
computed by integrating the cell density over the domain. If the cell density
function can be expressed analytically and the first order derivative of the
function exists, then a certain bin length $d$ is chosen and the cell count in
every bin of $d$ length is calculated. Then we generalize the center positions
of cells in every bin of length $d$, thus, the SPH approach can be
implemented, as it is indicated in Figure 4.2. However, it is not always
straightforward to obtain the analytical expression of cell density. If the
center positions of cells are given, the number of cells in each mesh element
can be counted, hence, the cell density will be computed analogously at each
mesh points. Therefore, the boundary value problem of cell density approach is
solved by numerical methods, for example, the finite-element methods.

ptCell densityCell count in length $d$Cell density approachCell positionsSPH
approach Figure 4.2: With exact expression of cell density function and the
first order derivative of the function exists, cell density approach is
implemented directly. Based on the cell density, the number of cells in a
certain region with length $d$ is determined and subsequently, the center
positions of cells can be generalized. Hence, the SPH approach is implemented.
ptCell positionsCell count in mesh elementSPH approachCell densityCell density
approach Figure 4.3: Given the center positions of cells, one can directly
implement the SPH model. Computing the number of cells in every mesh element
and divided by the length of the mesh element results into the cell density.
Subsequently, cell density approach can be implemented.

In this manuscript, all the numerical results are derived by finite-element
methods with Lagrangian linear basis functions. Regarding the first
implementation method (see Figure 4.2), we show the results with a Gaussian
distribution and sine function as cell density functions; see Figure 4.4 and
4.5. We start with the simulations in which we keep the number of cells and
the center positions of the cells the same, then we refine the mesh. In Figure
4.55(a)-5(c), the bin length $d$ is $0.35$, and the mesh size is a function of
$d$. The results solved by SPH approach become smoother. With various values
of $d$, the solutions to the approaches are overlapping only when the factor
between the $d$ and mesh size is closer to $1$. From Figure 4.55(d) to 5(f),
the mesh is fixed and we vary the value of $d$. We note that in Figure
4.55(f), the solution to the SPH approach is significantly different from the
solution to the cell density approach. It is mainly caused by the fact that
$d$ is too small and there is barely any fluctuation with the count of cells
in every $d$ length subdomain, while with the Gaussian distribution as the
cell density function, the majority of the cells are centered around $x=3.5$.
Hence, the solution to SPH approach still manages to be comparable with the
solution to the cell density approach; see Figure 4.44(f). Numerical results
of the simulation in Figure 4.4 are displayed in Table 4.1. There are some
noticeable differences between two approaches, in particular the convergence
rate in the $H^{1}$-norm: thanks to the given, differentiable cell density
function, the cell density approach converges faster. In addition, the cell
density approach requires less computational time with a factor of $15$.

(a) $h(d)=d/5,h=0.07,d=0.35$

(b) $h(d)=d/10,h=0.035,d=0.35$

(c) $h(d)=d/50,h=0.007,d=0.35$

(d) $d(h)=2h,h=0.07,d=0.14$

(e) $d(h)=h,h=0.07,d=0.07$

(f) $d(h)=h/4,h=0.07,d=0.0175$

Figure 4.4: The cell density function is Gaussian distribution and using the algorithm in Figure 4.2, different simulations are carried out with various mesh size and the total number of cells. Blue curves represent the solutions to $(BVP^{2}_{SPH})$, and ref curves are the solutions to $(BVP^{2}_{den})$ with $n_{c}(x)=50\times 1/\sqrt{2\pi\times 0.1^{2}}\exp\\{-(x-3.5)^{2}/(2\times 0.1^{2})\\}$. In Subfigure 4(a)–4(c), we set $d=0.35$ and cell positions are fixed, as $h$ is decreasing. From Subfigure 4(d) to Subfigure 4(f), we use the same finite-element method settings (where $h$ is sufficiently small with $h=0.07$), and simulations are carried out with various values of $d$. Table 4.1: Numerical results of two approaches in one dimension, where the cell density function is Gaussian distribution: $n_{c}(x)=50\times 1/\sqrt{2\pi\times 0.1^{2}}\exp\\{-(x-3.5)^{2}/(2\times 0.1^{2})\\}$. Here, we define $N_{s}=88$ and the mesh size $h=0.07$. The results are solved by finite-element method with algorithm in Figure 4.2.  |  SPH Approach |  Cell Density Approach
---|---|---
$\|u\|_{L^{2}((0,L))}$ |  $0.5441481069175041$ |  $0.36197930815501245$
$\|u\|_{H^{1}(((0,L))}$ |  $0.9642151731656272$ |  $0.871720645462775$
Convergence rate of $L^{2}-norm$ |  $1.75281178$ |  $1.826378221$
Convergence rate of $H^{1}-norm$ |  $0.66372231$ |  $1.816659924$
Reduction ratio of the subdomain $(a,b)$ (%) |  $13.88062$ |  $9.52381$
Time cost $(s)$ |  $0.045070$ |  $0.0032084$

(a) $h(d)=d/5,h=0.07,d=0.35$

(b) $h(d)=d/10,h=0.035,d=0.35$

(c) $h(d)=d/50,h=0.007,d=0.35$

(d) $d(h)=2h,h=0.07,d=0.14$

(e) $d(h)=h,h=0.07,d=0.07$

(f) $d(h)=h/4,h=0.07,d=0.0175$

Figure 4.5: The cell density function is sine function and using the algorithm
in Figure 4.2, different simulations are carried out with various mesh size
and the total number of cells. Blue curves represent the solutions to
$(BVP^{2}_{SPH})$, and ref curves are the solutions to $(BVP^{2}_{den})$ with
$n_{c}(x)=40|\sin(2x)|$. In Subfigure 5(a)–5(c), we set $d=0.35$ and cell
positions are fixed. From Subfigure 5(d) to Subfigure 5(f), we use the same
finite-element method settings (where $h$ is efficiently small with $h=0.07$),
and we take different values of $d$.

We consider cells that are located uniformly in the subdomain $(2,5)$, which
implies that the gradient or divergence of the cell density vanishes inside
the subdomain but does not exist at two endpoints of the subdomain. Hence, we
utilize the implementation method in Figure 4.3, as the center positions of
the cells are given, then the local cell density can be calculated per unit
area. Compared with the results shown in Figure 4.1, the results in Figure 4.6
and Figure 4.7 show the solutions to $(BVP^{2}_{SPH})$ and $(BVP^{2}_{den})$
respectively. Note that, in the finite-element method solutions, the magnitude
of the forces in both approaches are the same, and the variance of
$\delta_{\varepsilon}(x)$ is related to $h$ rather than $\Delta s$.
Furthermore, these figures verify that the convergence between SPH approach
and cell density approach is determined by the mesh size rather than by the
distance between any two adjacent cells. Table 4.2 displays the numerical
results of the simulation in Figure 4.6, in the perspective of the solution,
the reduction ratio of the subdomain and the computational cost. Similarly to
the figures, there is no significant difference between the norms and the
deformed length of the subdomain. However, the simulation time in the cell
density approach is much shorter than in the SPH approach with a factor of
$35$.

(a) $h=0.7$

(b) $h=0.07$

(c) $h=0.014$

(d) $h=0.007$

Figure 4.6: The finite-element method solutions to $(BVP^{2}_{SPH})$ and
$(BVP^{2}_{den})$ are shown where cells are uniformly located. With the fixed
positions of cells, the solutions are convergent as $h\rightarrow 0^{+}$. Blue
points are the centre positions of biological cells. Red curves represent the
solutions to $(BVP^{2}_{SPH})$ and blue curves represent the solutions to
$(BVP^{2}_{den})$.

(a) $N_{s}=10$

(b) $N_{s}=50$

(c) $N_{s}=100$

(d) $N_{s}=500$

Figure 4.7: The finite-element method solutions to $(BVP^{2}_{SPH})$ and $(BVP^{2}_{den})$ are shown with uniform distribution. Compared to the analytical result, the consistency between two approaches are unrelated to the number of cells, and the solutions are convergent as $h\rightarrow 0^{+}$. Here, we use $h=0.007$. Blue points are the centre positions of biological cells. Red curves represent the solutions to $(BVP^{2}_{SPH})$ and blue curves represent the solutions to $(BVP^{2}_{den})$. Table 4.2: Numerical results of two approaches in one dimension with biological cells located uniformly. Here, we define mesh size $h=0.07$ and $N_{s}=50$, which means $\Delta s=0.06$. The results are solved by finite-element method with algorithm in Figure 4.3. |  SPH Approach |  Cell Density Approach
---|---|---
$\|u\|_{L^{2}((0,L))}$ |  $0.1858655201$ |  $0.1858660118$
$\|u\|_{H^{1}((0,L))}$ |  $0.2780804415$ |  $0.2914497482$
Convergence rate of $L^{2}-norm$ |  $0.9940317098$ |  $0.9985295706$
Convergence rate of $H^{1}-norm$ |  $1.002001685$ |  $1.004380036$
Reduction ratio of the subdomain $(a,b)$ (%) |  $7.96908$ |  $7.98821$
Time cost $(s)$ |  $0.10391$ |  $0.0030458$

###  4.2 Two-Dimensional Results

In the multi-dimensional case, we are not able to write the analytical
solution to the boundary value problems. The results are all solved by the use
of the finite-element method applied to $(BVP^{3}_{SPH})$ and
$(BVP^{3}_{den})$. Note that the force magnitude of both boundary value
problems is the same. Following the same implementation methods as in one
dimension, simulations are carried out with two formulas of cell density: (1)
cells are located inside the subdomain $\Omega_{w}$ randomly by the uniform
distribution; (2) the cell density function is in the form of the standard
Gaussian distribution over the computational domain with
$n_{c}(\bm{x})=50\times\frac{1}{2\pi}\exp\left\\{-\frac{\|\bm{x}|^{2}}{2}\right\\}$.
Implementation methods in Figure 4.3 and Figure 4.2 are applied respectively
in Simulation (1) and (2).

For Simulation (1), no analytical expression for (the derivative of) the
density function is available. Therefore, the implementation starts with
generating the cell positions, according to the principles outlined in Figure
4.3. In Figure 4.8, the displacement results are displayed. From the figures,
hardly any significant differences between the solutions can be observed,
which indicates that these two approaches are numerically consistent. Table
4.3 displays more details about the two approaches regarding the numerical
analysis: most data are more or less the same. Thanks to the continuity of the
SPH approach, the convergence rate between two approaches is similar. However,
as it has been mentioned earlier, the agent-based model is computationally
more expensive than the continuum-based model; here, the difference is a
factor of $240$.

Figure 4.8: The displacement results are shown, which are solved from $(BVP^{3}_{SPH})$ and $(BVP^{3}_{den})$ when cells are randomly located in the subdomain $(-5,5)\times(-5,5)$. In other words, it is impossible to write the analytical expression of $n_{c}(\bm{x})$, subsequently, the algorithm in Figure 4.3 is selected. There are $196$ biological cells in the computational domain. Blue points are the center positions of biological cells, red curves are the original shapes of the subdomain, and the black curves represent the deformed boundary of the subdomain. Table 4.3: Numerical results of two approaches in two dimensions with random distribution for the positions of biological cells. Due to the nonexistence of divergence or gradient of cell density function, implementation method in Figure 4.3 is used. |  SPH Approach |  Cell Density Approach
---|---|---
$\|u\|_{L^{2}((0,L))}$ |  $10.105858093727422$ |  $12.314518769308366$
$\|u\|_{H^{1}((0,L))}$ |  $10.890035548437087$ |  $16.67100195050636$
Convergence rate of $L^{2}-norm$ |  $2.012144067$ |  $2.061050246$
Convergence rate of $H^{1}-norm$ |  $2.011995299$ |  $2.03760311$
Reduction ratio of the subdomain $\Omega_{w}$ (%) |  $24.24192$ |  $23.33667$
Time cost $(s)$ |  $4.20677$ |  $1.77219\times 10^{-2}$

According to the setting of the simulation, we define the cell density
function by

$n_{c}(\bm{x})=50\times\frac{1}{2\pi}\exp\left\\{-\frac{\|\bm{x}|^{2}}{2}\right\\},\mbox{in
$\Omega$},$

which is a Gaussian distribution multiplied by a positive constant. Similarly,
in two dimensions, the cell density is defined by the number of biological
cells per unit area. In other words, the cell count is computed by the local
cell density multiplied by the area of selected region. Here, we assume that
the selected region is a $1\times 1$ square, then we generate the center
positions of biological cells in every unit square based on the local number
of cells. Figure 4.9 shows the numerical results regarding two approaches.
There is no significant difference if the same mesh structure resolution is
used. As the mesh is refined, the solution to the SPH approach is smoother,
since the ”ring” in the center becomes more dominant. In Table 4.4, it can be
concluded that there are no significant differences between two approaches
except for the computational efficiency and the convergence rate of
$H^{1}$-norm. If the mesh is not fine enough, then the solution to the SPH
approach is less smooth, hence, the determination of the gradient of the
solution is less accurate.

(a) $h=0.6435$

(b) $h=0.3218$

(c) $h=0.1609$

Figure 4.9: The displacement results are shown, which are solved from $(BVP^{3}_{SPH})$ and $(BVP^{3}_{den})$ when cells are located according to the cell density function $n_{c}(\bm{x})=50\times\frac{1}{2\pi}\exp\\{-\frac{\|\bm{x}|^{2}}{2}\\}$. Hence, implementing algorithm in Figure 4.3 is used. There are $440$ biological cells in the computational domain. Blue points are the center positions of biological cells, red curves are the original shapes of the subdomain, and the black curves represent the deformed boundary of the subdomain. Table 4.4: Numerical results of two approaches in two dimensions with Gaussian distribution for the positions of biological cells. Figure 4.2 is implemented and there are $440$ biological cells in the computational domain. |  SPH Approach |  Cell Density Approach
---|---|---
$\|u\|_{L^{2}((0,L))}$ |  $11.14304909846569$ |  $13.188441094735877$
$\|u\|_{H^{1}((0,L))}$ |  $12.77795210802095$ |  $15.533928099123479$
Convergence rate of $L^{2}-norm$ |  $1.724918322$ |  $1.856132219$
Convergence rate of $H^{1}-norm$ |  $3.059831229$ |  $2.455660849$
Reduction ratio of the subdomain $\Omega_{w}$ (%) |  $19.65854$ |  $20.55949$
Time cost $(s)$ |  $0.62347$ |  $0.017315$

##  5 Conclusions

In this manuscript, we discussed the different models to simulate the pulling
forces exerted by the (myo)fibroblasts depending on different scales of the
wound region. We started from one dimension and later extended the models to
two dimensions. In one dimension, we can write explicitly the analytical
solution to the boundary value problem with specific distribution of the
locations of biological cells. The consistency of the solutions indicates the
possibility to prove the analytical connection between the two approaches. In
both one and two dimensions, the numerical solutions delivered by finite-
element methods with Lagrange linear basis functions implied that these two
models are consistent under certain mesh conditions (when the mesh size is
sufficiently small) and regardless the locations of the biological cells and
the implementation methods. In summary, regarding the displacement of the ECM
from the mechanical model, the agent-based model and the cell density model
are consistent from a computational point of view. This could be used to
transfer one type of model to the other one regarding the force balance in the
wound healing model, as the connection between these two models has been
suggested. We want to use the developed insights for the analysis of upscaling
between agent-based and continuum-based model formulations.

## References

  * Apostol and Ablow [1958] T. M. Apostol and C. Ablow.  Mathematical analysis.  _Physics Today_ , 11(7):32, 1958.
  * Boon et al. [2016] W. Boon, D. Koppenol, and F. Vermolen.  A multi-agent cell-based model for wound contraction.  _Journal of biomechanics_ , 49(8):1388–1401, 2016\.
  * Braess [2007] D. Braess.  _Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics_.  Cambridge University Press, 2007.
  * Cumming et al. [2009] B. D. Cumming, D. McElwain, and Z. Upton.  A mathematical model of wound healing and subsequent scarring.  _Journal of The Royal Society Interface_ , 7(42):19–34, 2009.
  * Enoch and Leaper [2008] S. Enoch and D. J. Leaper.  Basic science of wound healing.  _Surgery (Oxford)_ , 26(2):31–37, 2008.
  * Haberman [1983] R. Haberman.  _Elementary Applied Partial Differential Equations_ , volume 987.  Prentice Hall Englewood Cliffs, NJ, 1983.
  * Haertel et al. [2014] E. Haertel, S. Werner, and M. Schäfer.  Transcriptional regulation of wound inflammation.  In _Seminars in Immunology_ , volume 26, pages 321–328. Elsevier, 2014.
  * Koppenol [2017] D. Koppenol.  Biomedical implications from mathematical models for the simulation of dermal wound healing.  2017\.
  * Martin [1997] P. Martin.  Wound healing–aiming for perfect skin regeneration.  _Science_ , 276(5309):75–81, Apr. 1997.  doi: 10.1126/science.276.5309.75.  URL https://doi.org/10.1126/science.276.5309.75.
  * Olkin and Pratt [1958] I. Olkin and J. W. Pratt.  A multivariate tchebycheff inequality.  _The Annals of Mathematical Statistics_ , 29(1):226–234, Mar. 1958.  doi: 10.1214/aoms/1177706720.  URL https://doi.org/10.1214/aoms/1177706720.
  * Peng and Vermolen [2019a] Q. Peng and F. Vermolen.  Numerical methods to solve elasticity problems with point sources.  _Reports of the Delft Institute of Applied Mathematics, Delft University, the Netherlands_ , 1389-6520(19-02), 2019a.
  * Peng and Vermolen [2019b] Q. Peng and F. Vermolen.  Point forces and their alternatives in cell-based models for skin contraction.  _Reports of the Delft Institute of Applied Mathematics, Delft University, the Netherlands_ , 1389-6520(19-03), 2019b.
  * Peng and Vermolen [2020] Q. Peng and F. Vermolen.  Agent-based modelling and parameter sensitivity analysis with a finite-element method for skin contraction.  _Biomechanics and Modeling in Mechanobiology_ , 19(6):2525–2551, July 2020.  doi: 10.1007/s10237-020-01354-z.  URL https://doi.org/10.1007/s10237-020-01354-z.
  * Pukelsheim [1994] F. Pukelsheim.  The three sigma rule.  _The American Statistician_ , 48(2):88, May 1994\.  doi: 10.2307/2684253.  URL https://doi.org/10.2307/2684253.
  * Weisstein [2010] E. Weisstein.  Erf. from mathworld—a wolfram web resource.  _URL: http://mathworld.wolfram.com/Erf.html_ , 2010.