# Quality of approximating a mass-emitting object by a point source in a
diffusion model

Qiyao Peng1∗, Sander C. Hille1
1\. Mathematical Institute, Faculty of Science, Leiden University.
Neils Bohrweg 1, 2333 CA, Leiden, The Netherlands
* Correspondence<EMAIL_ADDRESS>
###### Abstract

For the sake of computational efficiency and for theoretical purposes, in
mathematical modelling, the Dirac Delta distributions are often utilized as a
replacement for cells or vesicles, since the size of cells or vesicles is much
smaller than the size of the surrounding tissues. Here, we consider the
scenario that the cell or the vesicle releases the diffusive compounds to the
immediate environment, which is modelled by the diffusion equation. Typically,
one separates the intracellular and extracellular environment and uses
homogeneous Neumann boundary condition for the cell boundary (so-called
spatial exclusion approach), while the point source approach neglects the
intracellular environment. We show that extra conditions are needed such that
the solutions to the two approaches are consistent. We prove a necessary and
sufficient condition for the consistency. Suggested by the numerical results,
we conclude that an initial condition in the form of Gaussian kernel in the
point source approach compensates for a time-delay discrepancy between the
solutions to the two approaches in the numerical solutions. Various approaches
determining optimal amplitude and variance of the Gaussian kernel have been
discussed.

Keywords: Diffusion equation, point source, Dirac delta distribution,
numerical analysis, Gaussian kernel

##  1 Introduction

In mathematical modelling, it is a common practice to replace a spatial object
with a negligible size by a point-particle, which has no volume. In classical
mechanics, objects may be replaced by point masses [16]. In electrostatics,
point charge [4] is a theoretical concept, used e.g. to describe the electric
field that results from a spatial distribution of charge. In biological
modelling, for instance, chemotaxis [13, 15] or in wound healing, when one
deals with a large-scale wound that is in the order of millimeter or even
centimeter, individual cells are regarded as particles and point sources that
secrete signalling molecules, which then diffuse in the surrounding
environment [7, 10].

These are all idealisations that are convenient not only for the modeller and
theorist. For computational efficiency of model simulation it may be even
necessary to employ point-particles. In particular when the objects are
moving, interacting or have internal dynamics of their own, operating at
several different scales. For example, simulation of several moving objects
that secrete a diffusing compound in a finite element approach can easily
become quite cumbersome [5, 12]. An implementation in terms of several moving
and ‘mass-emitting’ point-particles may then be better tractable [6], in
particular when the number of objects becomes large – but not so large that a
continuum description with densities or concentrations is a proper
representation.

What representation is proper, is determined by the underlying research
question for the modelling effort. It yields a tolerance for the deviation
from observations, which in turn yields a tolerance for deviations between
solutions of different modelling approaches. Moreover, in view of this
question and the implied tolerance it can be more or less relevant that a
point source cannot represent the spatial heterogeneity in shape or in flux
density of compound over the boundary of a truly spatial object.

Motivated by biological applications, in this paper we are concerned with the
mathematical question of the assessment of the quality of approximation when a
model for diffusion in the environment of a compound that is secreted by
several stationary spatial objects is replaced by that with point sources at
central locations of these objects. We do so by numerical simulations of two-
dimensional spatial configurations, in which the objects are equal, of
circular shape, with homogeneous constant flux density over their boundary
(see Figure 1.11(a) for the schematic set-up, for a single object). As such,
it can be viewed as a follow-up to an analytic approach in Evers et al. [3]
(for a single object) and a numerical approach in the setting of wound healing
dynamics in Peng and Vermolen [8, 9]. Heterogeneity in shape and flux density,
and movement of objects will be considered in follow-up research. If mass-
emitting point sources cannot approximate circular stationary objects with
homogeneous flux density over their boundary within tolerance, then this
cannot be expected either in the heterogeneous or in the non-stationary
setting.

(a) A schematic presentation of a single object ($\Omega_{C}$) within the
domain $\Omega$. $\Omega\setminus\Omega_{C}$ is the environment into which the
object emits a compound that freely diffuses there. It cannot escape through
the boundary $\partial\Omega$ of the domain. A typical direction of the normal
vector $\mathbf{n}$ on $\partial\Omega_{c}$ is indicated

(b) The concentration of compound in the environment of the spatial exclusion
approach at the last time step.

(c) The concentration of compound of the point source approach at the last
time step.

Figure 1.1: Schematic presentation (not in scale) of the spatial set-up in the
models, with circularly shaped stationary objects within a bounded domain
(Panel (a)). When there are multiple cells in the computational domain,
differences between the simulation results of the model with spatial objects
(Panel (b)) and that with mass-emitting point sources (Panel (c)) can be
observed. By proper choice of parameters and initial condition in the point-
source approach these can be kept within tolerance (see Section 3.3 for
simulation details). Here, the zero homogeneous initial condition is utilized,
which is pointing into the cell $\Omega_{C}$.

The paper aims at exhibiting in which cases and how – practically – a proper
representation can be made by means of mass-emitting point sources of the
spatial object setting. We propose a quantification of deviations between the
two approaches, find intolerable differences that occur and identify their
sources, and we propose how these can be resolved, e.g. concerning the
extension of initial conditions to a larger spatial domain (see Section 4 and
see Figure 1.11(b) and 1(c) for a visual comparison). We present a (mainly
numerical) analysis of consequences of particular choices for doing so on the
quality of approximation. We do not claim that the proposed choice of
extension by means of a Gaussian-shaped function is the best. It is a
reasonable and intuitively motivated choice for which we give an approach for
finding parameter settings for the extension that may keep the approximation
within accepted tolerance.

###  1.1 Mathematical formulation of the research questions

We consider a bounded domain $\Omega\subset\mathbb{R}^{2}$ with piece-wise
$C^{1}$-boundary $\partial\Omega$ in which there are embedded a finite number
$N_{c}$ of non-overlapping spatial objects, called cells, which are considered
as disjoint subdomains $\Omega_{C_{i}}$, ($i=1,\dots,N_{c}$), also with piece-
wise $C^{1}$-boundaries $\partial\Omega_{C_{i}}$, such that these boundaries
do not touch, nor touch the boundary $\partial\Omega$ of the initial domain.
Write $\Omega_{C}:=\bigcup_{i=1}^{N_{c}}\Omega_{C_{i}}$ for the totality of
cells. The complement $\Omega\setminus\overline{\Omega}_{C}$ will be called
the extracellular environment of the cells. For each cell we select a point
$\bm{x}_{c}^{i}\in\Omega_{C_{i}}$, which will function as center for
representing cell $C_{i}$ by a point-particle located at that point; see
Figure 1.11(a) for a schematic presentation of this set-up.

If $\bm{v}=(v_{1},v_{2})\in\mathbb{R}^{2}$, we write
$|\bm{v}|:=(|v_{1}|^{2}+|v_{2}|^{2})^{1/2}$ for its Euclidean norm.

####  1.1.1 Spatial exclusion model

The cells secrete a compound into the environment over their boundary with
prescribed flux density $\phi(\bm{x},t)\geq 0$ at
$\bm{x}\in\partial\Omega_{C}$ and time $t\geq 0$. This compound diffuses in
this environment according to Fickian diffusion with homogeneous diffusion
constant $D$, without further interaction. It cannot escape the domain
$\Omega$. Initially, there is a distribution $u_{0}(\bm{x})$ of this compound
in the environment. Thus, the density $u_{S}(\bm{x},t)$ at time $t$ of the
compound in the environment is described by the initial boundary value problem

$(BVP_{S})\quad\left\\{\begin{aligned} \frac{\partial
u_{S}(\bm{x},t)}{\partial t}-D\Delta u_{S}(\bm{x},t)&=0,&\mbox{in
$\Omega\backslash\bar{\Omega}_{C},t>0$,}\\\ D\nabla
u_{S}(\bm{x},t)\cdot\bm{n}&=\phi(\bm{x},t),&\mbox{on
$\partial\Omega_{C},t>0$,}\\\ D\nabla u_{S}(\bm{x},t)\cdot\bm{n}&=0,&\mbox{on
$\partial\Omega,t>0$,}\\\ u_{S}(\bm{x},0)&=u_{0}(\bm{x}),&\mbox{in
$\Omega\backslash\bar{\Omega}_{C}$,}\end{aligned}\right.$ (1.1)

where $\bm{n}$ is the outward pointing unit normal vector to the domain
boundary of $\Omega\setminus\bar{\Omega}_{C}$. Note that the flux density
$\phi(\bm{x},t)$ is positive at $\bm{x}\in\partial\Omega_{C}$ where there is
flux of compound into the environment $\Omega\setminus\bar{\Omega}_{C}$, as it
is shown in Figure 1(a).

As appropriate in the setting of the Finite Element Method (FEM) (cf. [14]) we
consider the weak solution concept for $(BVP_{S})$ in the spatial dimension,
which is given by

$(WF_{S})\left\\{\begin{aligned} &\text{Find $u_{S}(\bm{x},t)\in
H^{1}(\Omega\backslash\bar{\Omega}_{C})$, such that}\\\
&\int_{\Omega\backslash\bar{\Omega}_{C}}\frac{\partial
u_{S}(\bm{x},t)}{\partial
t}v_{1}(\bm{x},t)d\Omega+\int_{\Omega\backslash\bar{\Omega}_{C}}D\nabla
u_{S}(\bm{x},t)\cdot\nabla v_{1}(\bm{x},t)d\Omega\\\
&-\int_{\partial\Omega_{C}}\phi(\bm{x},t)v_{1}(\bm{x},t)d\Gamma=0,\\\
&\text{for any $v_{1}(\bm{x},t)\in
H^{1}(\Omega\backslash\bar{\Omega}_{C})$.}\end{aligned}\right.$

Here $d\Omega$ is the restriction of Lebesgue measure on $\mathbb{R}^{2}$ to
$\Omega$ and $d\Gamma$ denotes the surface measure on $\partial\Omega_{C}$,
that is so normalized that the Divergence Theorem holds without additional
constant.

####  1.1.2 Point source model

We want to compare the (weak) solution $u_{S}$ to $(BVP_{S})$ with the
solution to a suitable boundary value problem with point sources at the
locations $\bm{x}_{c}^{i}$ instead of spatial cells. These sources will be
expressed using Dirac measures $\delta_{\bm{x}_{c}^{i}}$ at $\bm{x}_{c}^{i}$,
or – equivalently – in the form of the Schwartzian delta distribution $\delta$
(see e.g. [11]), which is defined in any dimension by

$\langle\delta,f\rangle=f(\bm{0}),\qquad\mbox{if}\ f\in
C^{\infty}_{c}(\mathbb{R}^{n}).$

Informally written, as often done, the Dirac measure at $\bm{x}_{0}$ can then
be viewed as translation of the delta distribution:

$\int_{\Omega}f(\bm{x})d\delta_{\bm{x}_{0}}(\bm{x})\ =\
\int_{\Omega}\delta(\bm{x}-\bm{x}_{0})f(\bm{x})d\Omega=f(\bm{x}_{0}).$

The initial-boundary value problem defined by point sources is then given by

$(BVP_{P})\quad\left\\{\begin{aligned} \frac{\partial
u_{P}(\bm{x},t)}{\partial t}-D\Delta
u_{P}(\bm{x},t)&=\sum_{i=1}^{N_{c}}\Phi_{i}(t)\delta(\bm{x}-\bm{x}^{i}_{c}),&\mbox{in
$\Omega,t>0$,}\\\ D\nabla u_{P}\cdot\bm{n}&=0,&\mbox{on
$\partial\Omega,t>0$,}\\\ u_{P}(\bm{x},0)&=\bar{u}_{0}(\bm{x}),&\mbox{in
$\Omega,t=0$.}\end{aligned}\right.$ (1.2)

Here, $\Phi_{i}(t)$ is a function that describes the flux of mass per unit
time from the source at $\bm{x}_{c}^{i}$. We shall take

$\Phi_{i}(t)=\int_{\partial\Omega_{C_{i}}}\phi(\bm{x},t)d\Gamma.$ (1.3)

See Section 1.1.4 for further discussion of the selection of a suitable flux
function $\Phi_{i}$.

Again, we consider weak solutions in the setting of FEM. This amounts to the
following weak form – formulated for a single cell for convenience (dropping
indices $i$ that distinguish cells):

$(WF_{P})\left\\{\begin{aligned} &\text{Find $u_{P}(\bm{x},t)\in
H^{1}(\Omega)$, such that}\\\
&\int_{\Omega\backslash\bar{\Omega}_{C}}\frac{\partial
u_{P}(\bm{x},t)}{\partial
t}v_{2}(\bm{x},t)d\Omega+\int_{\Omega_{C}}\frac{\partial
u_{P}(\bm{x},t)}{\partial
t}v_{2}(\bm{x},t)d\Omega+\int_{\Omega\backslash\bar{\Omega}_{C}}D\nabla
u_{P}(\bm{x},t)\nabla v_{2}(\bm{x},t)d\Omega\\\ &+\int_{\Omega_{C}}D\nabla
u_{P}(\bm{x},t)\nabla
v_{2}(\bm{x},t)d\Omega=\int_{\Omega}\Phi(t)\delta(\bm{x}-\bm{x}_{c})v_{2}(\bm{x},t)d\Omega,\\\
&\text{for any $v_{2}(\bm{x},t)\in H^{1}(\Omega)$.}\end{aligned}\right.$

The singular nature of the delta distribution causes functional analytic
issues. Although weak solutions to $(WF_{F})$ exist and are unique (see Evers
et al. [3]), there does not exists a stationary solution in $H^{1}(\Omega)$.
The solution to the elliptic boundary value problem from the balance of
momentum with delta distribution is singular in the sense that for spatial
dimension higher than one, no formal solutions in the finite-element space
$H^{1}$ exist. Dealing with this singularity caused by delta distributions was
discussed in Peng and Vermolen [8, 9].

####  1.1.3 Measures for quantitative comparison

Quantifying the differences between the solutions to the spatial exclusion
model and the point source model in a meaningful way is not fully
straightforward. Various choices can be made. Hence, a first question is what
measure for comparison is meaningful in the context of applications. Here we
shall motivate those that we selected for use in this paper and one resulting
from Evers et al. [3].

First note that the two solutions ‘live’ on different spatial domains: $u_{P}$
is defined on $\Omega$, while $u_{S}$ is defined on the environment, the
subset $\Omega\setminus\bar{\Omega}_{C}$, only. Since there is no canonical
way of extending $u_{S}$ to the larger set $\Omega$, an objective comparison
of the two solutions is possible only on the smaller set
$\Omega\setminus\bar{\Omega}_{C}$. In accordance with this reasoning, the part
of solution $u_{P}$ on $\bar{\Omega}_{C}$ has not been shown in Figure
1.11(c).

The following proposition substantiates the intuition that any difference in
the two solutions is caused by a difference of flux over the boundary
$\partial\Omega_{C}$ of the two solutions, provided that their initial
conditions are the same on the environment of the cells.

######  Proposition 1.1.

Denote by $u_{S}(\bm{x},t)$ and $u_{P}(\bm{x},t)$ the weak solutions to the
spatial exclusion model $(BVP_{S})$ and the point source model $(BVP_{P})$,
respectively, and let $\partial\Omega_{C}$ be the boundary of the cells, from
which the compounds are released, with normal vector $\bm{n}$ pointing into
$\Omega_{C}$. Then

$\displaystyle\frac{1}{2}\frac{d}{dt}\bigl{\|}u_{S}-u_{P}\bigr{\|}^{2}_{L^{2}(\Omega\setminus\Omega_{C})}\
$ $\displaystyle=\
-D\int_{\Omega\setminus\Omega_{C}}\bigl{|}\nabla(u_{S}-u_{P})\bigr{|}^{2}d\Omega$
(1.4) $\displaystyle\qquad\qquad+\
\int_{\partial\Omega_{C}}(u_{s}-u_{P})(\phi-D\nabla
u_{P}\cdot\bm{n})\,d\Gamma.$

Assume moreover, that $u_{S}(\cdot,0)=u_{P}(\cdot,0)$ a.e. on
$\Omega\setminus\Omega_{C}$. Then, $u_{S}(\bm{x},t)=u_{P}(\bm{x},t)$ a.e. in
$\Omega\backslash\bar{\Omega}_{C}\times[0,\infty)$ if and only if

$\phi(\bm{x},t)-D\nabla u_{P}(\bm{x},t)\cdot\bm{n}=0,\qquad\mbox{ a.e. on
$\partial\Omega_{C}\times[0,\infty)$.}$

###### Proof.

See Appendix A. ∎

######  Remark 1.1.

The proposition can be analogously extended to any spatial dimension.

Therefore, the difference in total flux over the boundary, summed up to time
$t$, i.e.

$\int_{0}^{t}\int_{\Omega_{C}}\phi(\bm{x},s)-D\nabla
u_{P}(\bm{x},s)\cdot\bm{n}\;d\Gamma ds,$

yields a physically interpretable quantity of the deviation between the two
solutions. It is the difference in total amount of compound in the environment
between the spatial exclusion and point source approach, caused by the
difference in flux over the cell boundary $\partial\Omega_{C}$. Following the
definition of Evers et al. [3], we take the related quantity

$c^{*}(t):=\int_{0}^{t}\|\phi(\bm{x},s)-D\nabla
u_{P}(\bm{x},s)\cdot\bm{n}\|^{2}_{L^{2}(\partial\Omega_{C})}ds.$ (1.5)

as a measure of comparison. For technical reasons, we chose to work with
$L^{2}$-norm rather than the physically readily interpretable $L^{1}$-norm on
$\partial\Omega_{C}$. The former is easier accessible through the Finite
Element Method (FEM). Moreover, Evers et al. [3] gives various theoretical
estimates for $c^{*}(t)$, for a single cell and point source. Of course,
$\|f\|_{L^{1}(X,\mu)}\leq\mu(X)^{1/2}\|f\|_{L^{2}(X,\mu)}$.

The $L^{1}$-norm difference on the environment
$\Omega\setminus\bar{\Omega_{C}}$ compares the total amount of compound
between the two solutions. The $L^{1}$-norm difference of the gradients yields
information on differences in local fluxes that occur. For both we again
prefer to use the (related) $L^{2}$-norms, because of FEM that is used in the
numerical analysis.

We shall be looking for extension of the initial condition in the spatial
exclusion model to the cells $\Omega_{C}$, such that one arrives at an initial
condition on $\Omega$ for the point source approach that yields a good
approximation of the solution of the spatial exclusion model. To assess this
quality it is necessary to use a relative measure of comparison. That is, to
quantise a deviation in comparison to the total amount of compound in the
environment of the cells. A comparison of the effect of different initial
conditions with an absolute measure cannot be easily made, since different
conditions tend to give different amounts of compound in the environment.

####  1.1.4 Determining influx for the point source approach

When replacing a cell $C_{i}$ with mass flux density $\phi_{i}(\bm{x},t)$ over
its boundary $\partial\Omega_{C_{i}}$ in the direction of the environment by a
point source at $\bm{x}_{c}^{i}$ with mass flux $\Phi_{i}(t)$, one has to
decide how to relate the latter to the former. In Equation (1.3) we made the
choice that the total mass emitted by a cell at time $t$ per unit time is
equal to that emitted by the point source that replaces this cell. Thus, the
total mass emitted by the point source up to time $t$ is kept equal to the
total mass emitted by the cell into the environment. Since mass needs on
average a time of the order $\mathrm{diam}(C_{i})^{2}/4D$ to travel from the
center $\bm{x}_{c}^{i}$ to the boundary $\partial\Omega_{C_{i}}$, there will
be a time lag between the solution $u_{S}$ and $u_{P}$ on
$\partial\bar{\Omega}_{C}$.

Other choices for $\Phi_{i}(t)$ could be made. However, the definition in
Equation (1.3) seems most natural, particularly if one considers the
homogeneous flux. Moreover, the time lag between the solutions $u_{S}$ and
$u_{P}$ may be partially overcome by appropriately choosing the initial
condition $\bar{u}_{0}$ on $\Omega$ of solution $u_{P}$ in relation to the
initial condition $u_{0}$ on $\Omega\setminus\bar{\Omega}_{C}$ for $u_{S}$. In
this paper we shall focus on the latter means for minimizing the difference
between the two solutions.

####  1.1.5 Extension of initial condition

The major question addressed in this research is, in what way one can best
replace the initial condition $u_{0}$ for the spatial exclusion model by an
initial condition $\bar{u}_{0}$, defined on the whole domain $\Omega$, such
that the solutions $u_{S}$ and $u_{P}$ are ‘optimally close’. Here we fix the
flux relation as in Equation (1.3), as discussed in the previous section.

We consider $\bar{u}_{0}$ as an extension of $u_{0}$ to the larger domain. The
question is then, what function profile to take on $\Omega_{C}$, to make the
extension and how smooth does $u_{0}$ and this profile connects at the
boundary $\partial\Omega_{C}$?

In view of Proposition 1.1 and the further discussion in Section 1.1.3 it
seems reasonable to require that the flux over the boundary
$\partial\Omega_{C_{i}}$ in the point source model, created by the combined
effect of the selected initial condition profile in $\Omega_{C_{i}}$ and the
mass flux $\Phi_{i}(t)$ from the point source at $\bm{x}_{c}^{i}$ are as close
to $\phi_{i}(\bm{x},t)$ as possible over all time for which the solution is
computed. Again, various quantifiers for this difference may be selected.

###  1.2 Major assumptions

In this paper we limit our attention to a specific setting of the questions
raised above. First of all, we consider the general question of the quality of
approximating by a point source only in two spatial dimensions. This is mainly
to reduce computational intensity of the simulations. Moreover, the preceding
theoretical work [3] also considered a two-dimensional spatial domain.

The spatial objects are all taken circular in shape, with equal radius $r$. A
point source will emit compounds in a symmetric manner into its environment in
case of isotropic diffusion, which we consider here. One therefore knows a
priori, that replacing a non-circular object by a point source will not only
create error because of the reduction of the spatial object to a point, but
also because a non-symmetric object (for rotation) cannot be expected to
produce a symmetric emission profile, typically. The appropriate approximation
of a generally shaped object by one that is circular is another question, that
would be best considered separately.

The boundary $\partial\Omega_{C_{i}}$ of the circular cell $C_{i}$ of radius
$r$, centred at $\bm{x}_{c}^{i}$ is parameterized by the angle
$\theta\in[0,2\pi)$ relative to a reference direction. The parameterisation
$\gamma_{i}(\theta)$ is such that the measure $d\Gamma$ on
$\partial\Omega_{C_{i}}$, which is – recall – normalized such that the
Divergence Theorem holds without additional constants (see Section 1.1.1), is
given by

$\int_{\partial\Omega_{C_{i}}}f\,d\Gamma=r\int_{0}^{2\pi}f\bigl{(}\gamma_{i}(\theta)\bigr{)}\,d\theta,\qquad\mbox{for}\
f\in C(\partial\Omega_{C_{i}}).$ (1.6)

The simplest case of emission is one that is constant in time and spatially
homogeneous over the circular object. We consider that case here. That is,

$\phi(\bm{x},t)=\phi_{0}>0\qquad\mbox{for all}\ \bm{x}\in\partial\Omega_{C},\
t\geq 0.$ (1.7)

Spatially non-homogeneous flux over a circular boundary will allow to
approximate flux emitting from a non-circular object. Thus, this is certainly
an important aspect of the approximation question to consider. This setting is
currently being investigated. Results on this more general setting will appear
elsewhere.

We expect that the diffusion constant will influence the quality of
approximation. Therefore, we shall vary $D$ around a central value $D_{0}$ by
one order of magnitude.

###  1.3 Structure of the paper

The manuscript is structured as follows: The preparations for numerical
analysis and the approach taken are presented in Section 2. Furthermore, the
results of numerical simulations in Sections 3 and 4. The first considers
solutions with zero initial condition on the environment
$\Omega\setminus\Omega_{C}$, with non-zero extension to $\Omega$. The second
considers non-zero – but constant – initial condition on
$\Omega\setminus\Omega_{C}$ and suitable extension of this initial condition
to $\Omega$. Finally, conclusions and the outlook of this work are discussed
in Section 5.

##  2 Preparation for Numerical Analysis

For the numerical analysis we shall use non-dimensional versions of the models
$({BVP}_{S})$ and $(BVP_{P})$, in particular, their weak formulations.

###  2.1 Non-dimensional models

Denote the circular cell region centered at $\bm{x}_{c}$ and radius $r$ by
$\Omega_{C}:=\mathbb{B}(\bm{x}_{c},r)$. The entire domain is
$\Omega:=[-L,L]\times[-L,L]$. In the dimensionless model and simulations
thereof we scale space such that cell diameter becomes $1$. Thus, we get
spatial variables $\xi$ and computational domains given by

$\bm{\xi}:=\frac{\bm{x}}{2r},\
\hat{\Omega}_{C}:=\mathbb{B}\bigl{(}\bm{\xi}_{c},\displaystyle\frac{1}{2}\bigr{)},\
\hat{\Omega}:=\left[-\frac{L}{2r},\
\frac{L}{2r}\right]\times\left[-\frac{L}{2r},\ \frac{L}{2r}\right].$

Time is scaled by $\tau_{0}$ such that dimensionless time and the diffusion
constant become

$\tau:=\frac{t}{\tau_{0}},\quad\hat{D}:=\frac{D\tau_{0}}{4r^{2}}.$

Here, $\tau_{0}$ is chosen such that $\hat{D}=1$ corresponds to $D=D_{0}$, the
central value, and the relevant range for varying $D$ becomes
$\hat{D}\in[0.1,10]$. At last, we will consider $\phi(\bm{x},t)=\phi_{0}$
constant in time and space. We scale compound density by $u^{*}$ such that the
flux density in the new coordinates becomes $1$. That is,

$\frac{\phi_{0}\tau_{0}}{2ru^{*}}=1,\quad\gamma:=\frac{u}{u^{*}}.$

Then, $(BVP_{S})$ yields in a dimensionless system given by

$(BVP^{\prime}_{S})\left\\{\begin{aligned}
\frac{\partial\gamma_{S}(\bm{\xi},\tau)}{\partial\tau}-\hat{D}\Delta_{\bm{\xi}}\gamma_{S}(\bm{\xi},\tau)&=0,&\mbox{in
$\hat{\Omega}\backslash\bar{\hat{\Omega}}_{C},\tau>0$,}\\\
\hat{D}\nabla_{\bm{\xi}}\gamma_{S}(\bm{\xi},\tau)\cdot\bm{n}_{\bm{\xi}}&=1,&\mbox{on
$\partial\hat{\Omega}_{C},\tau>0$,}\\\
\hat{D}\nabla_{\bm{\xi}}\gamma_{S}(\bm{\xi},\tau)\cdot\bm{n}_{\bm{\xi}}&=0,&\mbox{on
$\partial\hat{\Omega},\tau>0$,}\\\
\gamma_{S}(\bm{\xi},0)&=\frac{u_{0}(\bm{\xi})}{u^{*}},&\mbox{in
$\hat{\Omega}\backslash\bar{\hat{\Omega}}_{C}$.}\end{aligned}\right.$

A similar transformation can be done analogously in the point source approach:

$(BVP^{\prime}_{P})\left\\{\begin{aligned}
\frac{\partial\gamma_{P}(\bm{\xi},\tau)}{\partial\tau}-\hat{D}\Delta_{\bm{\xi}}\gamma_{P}(\bm{\xi},\tau)&=\hat{\Phi}\delta(\bm{\xi}-\bm{\xi}_{c}),&\mbox{in
$\hat{\Omega},\tau>0$,}\\\
\hat{D}\nabla_{\bm{\xi}}\gamma_{P}\cdot\bm{n}_{\bm{\xi}}&=0,&\mbox{on
$\partial\hat{\Omega},\tau>0$,}\\\
\gamma_{P}(\bm{\xi},0)&=\frac{\bar{u}_{0}(\bm{\xi})}{u^{*}},&\mbox{in
$\hat{\Omega}$.}\end{aligned}\right.$

Here,

$\hat{\Phi}\ =\ \frac{\Phi\tau_{0}}{u^{*}\cdot 4r^{2}}\ =\ \frac{2\pi
r\phi_{0}\tau_{0}}{u^{*}\cdot 4r^{2}}\ =\ \pi,$ (2.1)

under the given scalings. The factor $4r^{2}$ in the denominator of the first
expression for $\hat{\Phi}$ in Equation (2.1) results from the different
behaviour under the scaling transformation $x\mapsto\xi=x/(2r)$ of Dirac
measure compared to Lebesgue measure in two dimensions, which is the reference
measure for the density functions on the left-hand side in the differential
equation.

Having made explicit the non-dimensional models ($BVP^{\prime}_{S}$) and
($BVP^{\prime}_{P}$) and having observed their similarity to the dimensional
models ($BVP_{S}$) and ($BVP_{P}$), we continue in our exposition with using
the original notation of $u_{S}$, $u_{P}$, etc. for the dimension-free
solutions and by omitting the ‘hats’ on parameters. Moreover, we stress that
in analysis and simulations we use the spatially homogeneous and constant
(rescaled) flux density, i.e. $\phi(\bm{x},t)=1$, in $(BVP_{S})$ and
corresponding flux at the point source $\Phi(t)=\Phi=\pi$, according to
Equation (2.1).

###  2.2 Computational Approach

In Evers et al. [3], the authors worked analytically on the upper bound of the
global difference between the solutions to these two approaches. Inspired by
it, a subsequent research question is ‘Under what circumstances and to what
extent, these two approaches are consistent?’ Therefore, the numerical
simulations are employed in this manuscript to quantify and visualize the
difference between the two solutions. A schematic representation of the mesh
structure of the two approaches are displayed in Figure 2.1. In other words,
for the spatial exclusion approach, we exclude the cell region from the
computational domain.

(a) The spatial exclusion approach

(b) The point source approach

Figure 2.1: A schematic representation of the two approaches and the mesh
structure. The computational domain is $(-10,10)\times(-10,10)$, and the cell
is located at $(-3.5,-4)$ with radius $\displaystyle\frac{1}{2}$. The blue
circles are the predefined cell boundary.

The parameter values that are used in this section are shown in Table 2.1, if
there is no further specification. Note that the parameters are dimensionless.
In this manuscript, finite-element methods and backward Euler are used for the
numerical simulations, for the spatial discretisation and time integration
respectively. Particularly, we use Python 3.10 and FEniCS package [1] version
2019.2.0.dev0. We bare in mind that in the implementation, instead of a smooth
circle, the cell region is constructed by a series of mesh points as a
polygon. In the point source approach, numerically, the initial condition at
the mesh points on the cell boundary $\partial\Omega_{C}$ is taken the same as
in the extracellular environment $\Omega\backslash\Omega_{C}$.

Table 2.1: Parameter values used in Section 3 and 4, corresponding to the dimensionless systems derived in Section 2.1 Parameter |  Value |  Description
---|---|---
$\hat{D}$ |  $0.1$ |  Diffusion coefficient
$L/(2r)$ |  10 |  Size of the computational domain
$\Delta\tau$ |  $0.04$ |  Time step
$T$ |  $40$ |  Total time
$h$ |  $0.127$ |  Average mesh size

###  2.3 Preliminary Results: Effects of Varying the Diffusion Constant

A diffusing particle that is released at the centre of a circular cell of
radius $R=1/2$ reaches the cell boundary on average on a time scale $R^{2}/D$.
When the time step $\Delta\tau$ is larger than $R^{2}/D$, the global
difference between the norms of the solutions in both approaches is not really
significant; see Figure 2.22(a)-2(b). This is mainly due to the fact that
within one time step, the compounds have already reached the boundary of the
cell in the point source approach. Subsequently, the diffusion basically
starts from the boundary of the cell at $t=\Delta\tau$. However, when we
decrease the diffusion coefficient significantly, longer time is needed for
the compounds to reach the boundary of the cell in the point source approach.
In other words, the influx from the point source firstly needs to ‘fill’ the
intracellular space $\Omega_{C}$, before it can mimic diffusion from the
boundary of the cell.

(a) $D=10$

(b) $D=10$

(c) $D=1$

(d) $D=1$

(e) $D=0.1$

(f) $D=0.1$

Figure 2.2: Several simulations were conducted with varying diffusion
coefficients, namely, $D\in\\{10,1,0.1\\}$. Zero initial conditions were taken
in all cases. In the subfigures in the left column, various global norms (i.e.
$\|u\|_{L^{2}(\Omega\backslash\Omega_{C})},\|u\|_{H^{1}(\Omega\backslash\Omega_{C})},\|\nabla
u\|_{L^{2}(\Omega\backslash\Omega_{C})}$) are shown for both approaches. In
the subfigures in the right column, we show the local norm differences of all
the categories of aforementioned norms. The increase in $L^{2}$-norm of the
solution on $\Omega\setminus\Omega_{C}$ – hence also in $H^{1}$-norm – in (a),
(c) and (e) is due to the constant influx of compound in both spatial
exclusion and point source approach. Here, $H^{1}(\Omega\setminus\Omega_{C})$
is equiped with the Hilbert space norm
$\|f\|_{H^{1}}^{2}=\|f\|_{L^{2}}^{2}+\|\,|\nabla f|\,\|_{L^{2}}^{2}$.

Notice that the local difference between the two solutions increases in all
the norms with decreasing diffusivity $D$; see Figure 2.2 2(b), 2(d) and 2(f).
This is caused by the effect that even though compounds may have reached the
boundary of the cell in the point source approach after some time, these can
never fully ‘catch up’ with the amount of mass already released by this
boundary from the start in the spatial exclusion approach. Thus, one gets a
systematic delay, visible in Figure 2.2.

Furthermore, as time proceeds, all the norms of differences slowly decrease
for $D=1$ and $D=0.1$. The $L^{2}$-norm is controlled by expression in
Equation (1.4). After a transient, the difference in flux density over the
boundary $\partial\Omega_{C}$ becomes small; see e.g. Figure 1(c) that
indicates that concentration differences are small, hence also the differences
in flux. Thus, the boundary-integral term in Equation (1.4) converges to 0.
The first term, with the difference in gradients of solutions, will also
converge to 0. It depends on the speed of convergence of both terms, whether
the sign of the derivative of the $L^{2}$-norm of the difference will be
positive or negative. This phenomenon can be observed in Figure 2.2 2(b), for
$D=10$, where the $L^{2}$-norm of the difference is increasing towards the end
of the simulation time interval.

This phenomenon seems to persist with decreasing $D$ until it disappears for
$D$ between $1$ and $5$. To check whether this is a phenomenon in the
numerical simulation rather than the model, we decreased the time step
substantially. This did not affect the computed $L^{2}$-norm. We conclude,
that it is a consequence of a change in dominance between the two integral
terms in Equation (1.4). It is not yet clear what controls this change
precisely.

For the remainder of this paper, we shall focus on a second observation that
can be made from the preliminary simulations: varying $D$, these show that for
smaller diffusion constant, there is larger ‘delay’ between the solutions, see
Figure 2(e). This causes a systematic deviation between the two solutions.
Hence, it is important to see how to reduce this effect. Therefore, in all
following numerical simulations, we select $D=0.1$, since it gives the largest
deviation in norms between the spatial exclusion model and the point source
model.

##  3 Extending Zero environmental Initial Value

We shall first examine the issue of systematic delay for the point source
solution when there is a zero initial condition in the environment for the
spatial exclusion model. The situation where the initial condition there is
constant, but non-zero, turned out to require a modified approach. It is
discussed in the next section.

Since the solutions in the two approaches are defined on different spatial
domains, the initial condition for the spatial exclusion model cannot be
simply carried over to the larger domain of the point source model: choices
for extension to the complement $\bar{\Omega}_{C}$ must be made. In the
preliminary results of Section 2.3 we have seen that simply extending by zero
leads to structural delay in the solution, especially for small diffusion.

###  3.1 A Gaussian-shaped extension

We consider now a single cell, centered at $\bm{x}_{c}$. According to
Proposition 1.1, it is expected that the initial condition in $(BVP_{P})$,
i.e. $\bar{u}_{0}(\bm{x})$, cannot be simply set to zero on $\Omega_{C}$ as
extension of the zero initial value $u_{0}(\bm{x})$ of $(BVP_{S})$ on
$\Omega\setminus\Omega_{C}$. As the spatial exclusion approach removes the
cell region in the computational domain and starts diffusion from the cell
boundary directly, while in the point source approach, firstly the compound
needs to reach the cell boundary from the cell center, which takes extra time.
To compensate this time difference and inspired by the fact that the
difference between two approaches (see Figure 2.22(c)-2(f)) stays more or less
constant, we set $\bar{u}_{0}(\bm{x})$ on $\Omega_{C}$ in the form of the
fundamental solution to the diffusion equation on $\mathbb{R}^{d}$, which is
given by [2]

$P^{D}(\bm{x},t)=\left\\{\begin{aligned} &\frac{1}{(4\pi
Dt)^{d/2}}\exp\left\\{-\frac{|\bm{x}|^{2}}{4Dt}\right\\},&\mbox{
$t>0,\bm{x}\in\mathbb{R}^{d}$,}\\\ &0,&\mbox{
$t<0,\bm{x}\in\mathbb{R}^{d}$,}\end{aligned}\right.$ (3.1)

where $d$ is the dimension. In this study, we consider $d=2$ only. If we
imagine that diffusion inside the cell has started a time $t_{0}>0$ before the
start of the point source model from a unit Dirac mass at the centre
$\bm{x}_{c}$, then at time $t$ it will have reached a distribution on
$\mathbb{R}^{2}$ of Gaussian shape

$P_{t+t_{0}}^{D}(\bm{x},\bm{x}_{c})=\frac{1}{4\pi
D(t+t_{0})}\exp\left\\{-\frac{|\bm{x}-\bm{x}_{c}|^{2}}{4D(t+t_{0})}\right\\},\qquad\mbox{
$t,t_{0}>0,\ \bm{x}\in\mathbb{R}^{2}$.}$ (3.2)

We can modify $t_{0}$ and the intensity $p_{0}>0$ of the initial condition to
arrive at a proposed extension of the initial condition $\bar{u}_{0}(\bm{x})$
in $(BVP_{P})$ as a (discontinuous) truncated and scaled fundamental solution:

$\bar{u}_{0}(\bm{x})=\begin{cases}\
p_{0}P_{t_{0}}^{D}(\bm{x},\bm{x}_{c}),\quad&\bm{x}\in\bar{\Omega}_{C},\\\ \
0,&\bm{x}\in\Omega\backslash\bar{\Omega}_{C}.\end{cases}$ (3.3)

The idea is, to choose $(p_{0},t_{0})$ in such a way that the flux condition
of Proposition 1.1 is met in the best possible way.

Denoting $r=|\bm{x}-\bm{x}_{c}|$, the distance to the singular point, the flux
density at $\bm{x}$ in the direction pointing away from $\bm{x}_{c}$ that
originates from the initial condition only is given by

$\phi_{1}(r,t)=-D\nabla(p_{0}P_{t+t_{0}}^{D}(r))\cdot\bm{n}=-Dp_{0}\frac{\partial
P_{t+t_{0}}^{D}(r)}{\partial r}=\frac{p_{0}r}{2(t+t_{0})}P_{t+t_{0}}^{D}(r).$

Hence, the flux density over $\partial\Omega_{C}$ from the initial condition
in Equation (3.3) reads as

$\phi_{1}(R,t)\ =\ \frac{p_{0}R}{2(t+t_{0})}P_{t+t_{0}}^{D}(R)\ =\
\frac{p_{0}R}{8\pi
D(t+t_{0})^{2}}\exp\left\\{-\frac{R^{2}}{4D(t+t_{0})}\right\\}.$ (3.4)

With the production of compounds at the center of the cell $\bm{x}_{c}$ and
production rate $\Phi(\bm{x}_{c})$, the fundamental solution of $(BVP_{P})$ is
given by

$\displaystyle u_{P}(\bm{x},t)$
$\displaystyle=\int_{0}^{t}\int_{\mathbb{R}^{2}}P^{D}_{t-s}(\bm{x},\bm{y})\delta(\bm{y}-\bm{x}_{c})\Phi(\bm{x}_{c})d\bm{y}ds$
$\displaystyle=\int_{0}^{t}\Phi(\bm{x}_{c})P^{D}_{t-s}(\bm{x},\bm{x}_{c})ds\
=\ \int_{0}^{t}\Phi(\bm{x}_{c})P^{D}_{t-s}(r)ds.$

Subsequently, the flux density caused by the point source at $\bm{x}_{c}$ only
is computed as

$\displaystyle\phi_{2}(r,t)$ $\displaystyle=-D\nabla u(r,t)\cdot\bm{n}\ =\
-D\frac{\partial u(r,t)}{\partial r}$ $\displaystyle=D\frac{\partial}{\partial
r}\int_{0}^{t}\Phi(\bm{x}_{c})P^{D}_{t-s}(r)ds\ =\
\frac{\Phi(\bm{x}_{c})}{2\pi r}\exp\left\\{-\frac{r^{2}}{4Dt}\right\\}.$

Then, we obtain the flux over $\partial\Omega_{C}$ as

$\phi_{2}(R,t)=\frac{\Phi(\bm{x}_{c})}{2\pi
R}\exp\left\\{-\frac{R^{2}}{4Dt}\right\\}.$ (3.5)

According to Proposition 1.1, and given $-D\nabla
u_{P}\cdot\bm{n}=\phi_{1}(R,t)+\phi_{2}(R,t)$ over $\partial\Omega_{C}$, the
relation between the scale $p_{0}$ and prediffused time $t_{0}$ can be
determined from the approximate equation

$\displaystyle\phi(\bm{x},t)\ $ $\displaystyle=\ -D\nabla u_{P}\cdot\bm{n}\
\approx\ \phi_{sum}:=\ \phi_{1}(R,t)+\phi_{2}(R,t)$ $\displaystyle\ =\
\frac{p_{0}R}{8\pi
D(t+t_{0})^{2}}\exp\left\\{-\frac{R^{2}}{4D(t+t_{0})}\right\\}\ +\
\frac{\Phi(\bm{x}_{c})}{2\pi R}\exp\left\\{-\frac{R^{2}}{4Dt}\right\\}.$ (3.6)

Let $t=0$, then

$p_{0}(t_{0})\approx\frac{2t_{0}\phi(\bm{x},t)}{RP^{D}_{t_{0}}(R)}.$ (3.7)

###  3.2 A Comparison: Gaussian-shaped Initial Value inside the Cell

We set $D=0.1$, such that the time delay between the solutions in the two
approaches is the largest in the range of $D$ that we consider, see Section
2.3. We take as initial condition $\bar{u}_{0}$ on $\Omega$ for the point
source model the function defined in Equation (3.3). Strictly speaking, this
function is not in $H^{1}(\Omega)$. However, numerically the function value of
$\bar{u}_{0}$ at the mesh points in the interior of the cell, i.e. in
$\Omega_{C}$ (and 0 for those on the boundary $\partial\Omega_{C}$, see
Section 2.2) are used to obtain a numerical approximation for $\bar{u}_{0}$
through FEM that is in $H^{1}$. This is the extension $\bar{u}_{0}$ that is
actually considered, but which is difficult to prescribe explicitly.

Thus, one has two degrees of freedom in selecting the Gaussian-shaped initial
condition. The question is, how to choose the value of $(p_{0},t_{0})$. As
first option we considered determining $(p_{0},t_{0})$ by minimizing the total
deviation between $\phi_{sum}(R,t)$ and $\phi(\bm{x},t)=1$ over the time
interval $[0,T]$ in $L^{1}$-sense. That is, $(p_{0},t_{0})$ is taken as

$(\hat{p}_{0},\hat{t}_{0})\in\operatorname*{argmin}_{(p_{0},t_{0})}\int_{0}^{T}|\\\
Phi_{sum}(t)-\phi(\bm{x},t)|dt,$ (3.8)

This choice is referred to as ‘Option 1’.

Figure 3.1 shows the global norms of the solutions and norms of differences
between the spatial exclusion approach and the point source approach with two
types of initial condition: zero initial condition compared to a Gaussian-
shaped initial condition on the cell’s interior. The shape parameter
$(p_{0},t_{0})$ has been determined according to Option 1, i.e. Equation
(3.8). In Figure 3.1(a) one sees that the norms of the individual solutions of
spatial exclusion and point source model will converge to each other after a
transient, when the Gaussian-shaped extension of the initial condition is
used. It indicates that the use of Gaussian-shaped initial condition does
compensate the time delay, compared to the use of zero initial condition. The
norms of the difference of the two solutions (with Gaussian-shaped extension)
converge quickly to a steady – though non-zero – level. So, error is
controlled well after an initial transient. Moreover, the error is less than
that in the case of the zero-extension, as can be seen in Figure 3.1(b). Thus,
the Gaussian-shaped extended initial condition yields an improvement of the
approximation over the homogeneously-zero extended initial condition.

When the Gaussian-shaped initial condition is used, due to the discontinuity
of $\bar{u}_{0}(\bm{x})$ given by Equation (3.3) in the point source approach,
the $H^{1}$-approximation of $\bar{u}_{0}(\bm{x})$ that is taken as initial
condition in the numerical analysis has a large gradient locally in a small
region of the boundary of the cell. This produces a spike of
$\|u_{S}-u_{P}\|_{H^{1}(\Omega\backslash\Omega_{C})}$ and $\|\nabla
u_{S}-u_{P}\|_{L^{2}(\Omega\backslash\Omega_{C})}$ that appear in the
beginning of the simulation. However, the differences quickly dropped to a
lower level, compared to when the homogeneously-zero initial condition is
used.

(a) Global norm differences of solutions to the two approaches

(b) Local norm differences of solutions to the two approaches

(c) $c^{*}(t)$

Figure 3.1: For $D=0.1$, computed (a) global norms for solutions and (b) local
norm differences (namely, the $L^{2}-$ and $H^{1}-$norm of the solution $u$,
and the $L^{2}-$norm of the gradient of the solution $u$) between the two
approaches, as well as (c) the quantity $c^{*}(t)$. Here, the amplitude and
the variance of the Gaussian-shaped initial condition are determined by
minimizing the total flux deviation on $\partial\Omega_{C}$ over $[0,T]$, i.e.
$(p_{0},t_{0})$ have been taken as in Option 1 in Table 3.1. Black curves
represent results from the spatial exclusion approach, and blue and red curves
represent the zero and Gaussian-shaped initial conditions, respectively.
Solid, dashed and dotted lines represent
$\|u_{S}-u_{P}\|_{L^{2}(\Omega\backslash\Omega_{C})}$,
$\|u_{S}-u_{P}\|_{H^{1}(\Omega\backslash\Omega_{C})}$ and $\|\nabla
u_{S}-\nabla u_{P}\|_{L^{2}(\Omega\backslash\Omega_{C})}$, respectively. Note
that the introduction of the Gaussian-shaped initial condition in the point
source approach removed the time delay, essentially.

Figure 3.11(c) shows the time-integrated deviation between the prescribed flux
$\phi(\bm{x},t)=1$ in $(BVP_{S})$ and the flux generated from $(BVP_{P})$ on
the cell boundary $\partial\Omega_{C}$ (i.e. Equation (1.5)). Due to the
discontinuity in the initial condition for $(BVP_{P})$, the gradient of the
flux in $(BVP_{P})$ over $\partial\Omega_{C}$ is large, hence, the $c^{*}(t)$
is larger when the Gaussian-shaped initial condition is used than for the zero
initial condition. As shown in Figure 3.11(c), this measure of quality reaches
a rather small steady rate of increase quicker for the Gaussian-shaped
condition. However, upon quick inspection, 3.11(c) seems to indicate that the
zero initial condition performs better than the Gaussian-shaped extension,
which is not the case in view of the norms of local differences presented in
Figure 3.1(b).

We conclude that the Gaussian-shaped extension of the initial condition
improves the quality of approximation compared to the homogeneously-zero
initial condition. However, the graph of $c^{*}(t)$ turns out to be difficult
to interpret towards drawing conclusions on the quality of approximation in
terms of $L^{2}$\- and $H^{1}$-norms.

###  3.3 Multiple Cells in the Computational Domain

The chief benefits of the point source model over the spatial exclusion model
lie in analytical tractability and computational efficiency when there are
many cells, in particular when these are also moving. Therefore, we shall now
examine effects that may occur due to the presence of multiple cells. That
there will be some effect of multiple cells can be anticipated. Intuitively
speaking, part of the flux of mass released by one cell in the point source
model will at some point in time freely transverse the part of the domain that
is the interior of another cell in the spatial exclusion model. In the latter,
diffusing particles would have reflected on the boundary of this cell. So,
different trajectories of the diffusion process are expected, which may result
in differences in the solutions. The impact of this phenomenon will depend on
the distance between the cells, their size and the diffusivity.

In this study, we shall only exhibit the impact of this phenomenon.
Investigation of the question how to compensate this impact is deferred to
another study. For simplicity, we consider two and ten cells respectively with
different locations. All the cells are assumed to be identical in shape, size
and (constant) flux density over the boundary. The locations of the cells over
the domain are shown in Figure C.1 in Appendix C. Similarly to the previous
section, we plot against time various norm differences of the solutions, and
$c^{*}(t)$ of the same single cell, i.e. the cell that is present at the same
location in all configurations. In every subfigure in Figure 3.2, the left
panel shows the result when the homogeneous initial condition is applied, i.e.
$u_{0}(\bm{x},0)=0$ for the computational domain, and the right panel is when
the Gaussian-shaped extension is used as the initial condition inside the cell
$\Omega_{C}$.

(a) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(b) $\|\nabla u_{S}(\bm{x},t)-\nabla
u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(c) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{H^{1}(\Omega\backslash\Omega_{C})}$

(d) $c^{*}(t)$ of the same cell

Figure 3.2: Numerical results of multiple cells in the computational domain.
Cells are identical ($\phi_{i}=1$, $R_{i}=\frac{1}{2}$). $D=01.$. Results have
been plotted for homogeneously-zero initial condition (solid curve in the
left) and Gaussian-shaped initial condition (dashed curve in the left) in
every subfigure. Different colors of curves represent different number of
cells, as indicated. The used cell configurations for the computations are
shown in Figure C.1.

Generally speaking, we observe the existence of the impact from other cells,
since for every error quantifier, the more cells in the computational domain,
the larger the quantity. Furthermore, the distance between the cells plays a
significant role in this impact: for the two-cell case, when the cells are
closer to each other, the influence appears sooner, and $c^{*}(t)$ is
increasing continuously and seems to not be able to reach a steady state,
which is opposite to the two-cell case with a larger distance. As the number
of cells increases, the quantities are increasing as well except for the
$c^{*}(t)$, which results from the minimal distance between the cells being
smaller than the two-cell case with a smaller distance. In other words, it
verifies the assumption that the distance between cells is a significant
factor to study the consistency of the solutions to the two approaches when
there are multiple cells in the computational domain. Also in this setting, a
Gaussian-shaped extension of the initial conditions yields a better
approximation, although the difference tends to become smaller when the number
of cells increases.

###  3.4 Optimal Selection of $(p_{0},t_{0})$ in the Initial Condition

So far, we compute the amplitude $p_{0}$ and variance parameter $t_{0}$ by
minimizing expression (3.8). However, there are other options to select
$(p_{0},t_{0})$ by minimizing different objective functions. Key objective is
to minimize the deviation between $\phi(\bm{x},t)$ and $\phi_{sum}(t)$ as
expressed in Equation (3.6). The $L^{1}$-norm deviation over time, as in
Equation (3.8) may be replaced by e.g. supremum norm. In Lemma B.1 we
summarized characteristic properties of $\phi_{sum}(t)$. It may have at most
two extreme values, one above the (constant) level $\phi$, one below. In order
to make both extreme deviations from the target value $\phi$ as small as
possible, one can take as objective to minimize

$\bigl{|}\max_{0\leqslant t\leqslant T}(\phi_{sum}(t)-\phi)\bigr{|}\ +\
\bigl{|}\min_{0\leqslant t\leqslant T}(\phi_{sum}(t)-\phi)\bigr{|},$ (3.9)

over $(p_{0},t_{0})$. One may think too of starting the approximation
$\phi_{sum}$ at level $\phi$ at $t=0$. This yields Equation (3.7) as
constraint on the $(p_{0},t_{0})$ value pairs. There is no a priori guarantee
however, that this yields the best result for the corresponding unconstrained
minimisation problem with the same objective function. A combination of
objective functions defined in Equations (3.8) and (3.9) has also been
considered.

Table 3.1 gives an overview of the objective functions that have been
examined, each with and without the value-pair constraint (Equation (3.7)).
The $(p_{0},t_{0})$-values thus obtained for each option, by applying the
optimize.minimize function of Python package Scipy (version 1.8.0) are shown
in the last two columns.

Table 3.1: Considered options to compute $(p_{0},t_{0})$ as minimisation of the indicated different objective functions and constraint on the $(p_{0},t_{0})$-value pair. Initial condition is $u_{0}=0$ in $\Omega\backslash\Omega_{C}$. The optimization was performed by the function optimize.minimize of Scipy package (version $1.8.0$) in Python. Options |  Objective function |  Constraints |  Value of $p_{0}$ |  Value of $t_{0}$
---|---|---|---|---
Option 1 |  $\int_{0}^{T}|\phi_{sum}(t)-\phi|dt$ |  - |  $55.379$ |  $3.660$
Option 2 |  |  $\displaystyle p_{0}(t_{0})=\frac{2t_{0}\phi}{RP^{D}_{t_{0}}(R)}$ |  $37.114$ |  $2.383$
Option 3 |  $|\max_{0\leqslant t\leqslant T}(\phi_{sum}(t)-\phi)|$ |  - |  $27.946$ |  $1.905$
Option 4 |  $+|\min_{0\leqslant t\leqslant T}(\phi_{sum}(t)-\phi)|$ |  $\displaystyle p_{0}(t_{0})=\frac{2t_{0}\phi}{RP^{D}_{t_{0}}(R)}$ |  $21.737$ |  $1.737$
Option 5 |  $\int_{0}^{T}|\phi_{sum}(t)-\phi|dt+|\max_{0\leqslant t\leqslant T}(\phi_{sum}(t)-\phi|$ |  - |  $31.451$ |  $2.086$
Option 6 |  $+|\min_{0\leqslant t\leqslant T}(\phi_{sum}(t)-\phi)|$ |  $\displaystyle p_{0}(t_{0})=\frac{2t_{0}\phi}{RP^{D}_{t_{0}}(R)}$ |  $34.439$ |  $2.283$

Figure 3(d) shows the plot of the important quantifiers of the consistency
between the approaches. As a reference, we add the curve which is obtained by
defining zero initial conditions over $\Omega$ in $(BVP_{P})$. Extending the
initial condition on $\Omega_{C}$ with the Gaussian kernel provides sufficient
flux to reduce the local difference between the solutions in these two
approaches, which can be seen in Figure 3(d)(a). Nevertheless, due to
discontinuity of the initial condition on $\partial\Omega_{C}$ as defined in
Equation (3.3), similarly to Figure 3.1(a), using the inhomogeneous initial
condition causes that the $H^{1}-$norm (see Figure 3(d)(d)) and the
$L^{2}-$norm of the gradient of the solution (see Figure 3(d)(c)) difference
are significantly larger than using the homogeneous initial condition in the
first few time steps. However, when these have stabilized after this transient
time period, the deviations for the all Gaussian-shaped extensions of the
initial condition are almost one order of magnitude lower than those for the
homogeneously-zero initial condition.

Again, the quantifier $c^{*}(t)$, defined in Equation (1.5), does not allow
for a conclusion as clear as that provided by the $L^{2}$ and $H^{1}$-norm
differences. It also contains the gradient of $u_{P}$. Hence, most graphs of
$c^{*}(t)$ in Figure 3(d)(d) are above the graph for the homogeneously-zero
initial condition (dashed curve), except for Option 3 and Option 4. There, the
maximal and minimal difference between the presumed flux density in
$(BVP_{S})$ and the analytical flux $D\nabla u_{P}\cdot\bm{n}$ from
$(BVP_{P})$ are minimized. In other words, Option 3 and 4 select
$(p_{0},t_{0})$ such that mostly $D\nabla u_{P}\cdot\bm{n}$ is close to
$\phi(\bm{x},t)$. In particular, Option 4 appears to be the best option, since
$c^{*}(t)$ reaches the steady state fastest and with the smallest value. In
the latter option one enforces that $\phi(\bm{x},0)=\phi_{sum}(0)$. However,
Figure 3(d)(a) indicates that Option 4 is the worst choice among the Gaussian-
shaped extensions, from the point of view of the $L^{2}$-norm distance. There
is no clear distinction among the options from the point of view of the
$H^{1}$-norm.

(a) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(b) $\|\nabla u_{S}(\bm{x},t)-\nabla
u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(c) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{H^{1}(\Omega\backslash\Omega_{C})}$

(d) $c^{*}(t)$

Figure 3.3: Measures of quality of approximation between solutions as a
function of time. The Gaussian-shape parameters $(p_{0},t_{0})$ have been
computed according to the options listed in Table 3.1 and the resulting values
for $(p_{0},t_{0})$ have been used in the simulation
($D=0.1,\phi=1,R=\frac{1}{2}$). Blue dashed line represents the $(BVP_{P})$
when $u_{P}(\bm{x},0)=0$ for $\Omega$, and other coloured solid lines
represent the $(BVP_{P})$ when the initial condition is given by Equation
(3.3). The curve for Option overlaps visually with that of Option 2 in Panel
(a), (b) and (c)

.

In view of Proposition 1.1, Equation (1.4), the behaviour of the $L^{2}$ norm
is a delicate interplay between the flux difference over the boundary
(measured by $c^{*}(t)$ and the difference of the gradient of the solutions.
These two quantities are not independent. More analytic insight – if
obtainable – is required to understand properly the apparent discrepancy
between the $c^{*}$ measure and $L^{2}$-norm of difference of solutions.

In general, the smaller the quantifier is, the better the option. However,
again there seems to be no “best” option for all quantifiers considered.
Again, it turns out that $c^{*}(t)$ has to be interpreted with greatest care.

##  4 Extending Nonzero environmental Initial Value

As observed in the previous results, discontinuity over $\partial\Omega_{C}$
caused by the Gaussian-shaped initial condition in the interior of the cell
and the zero homogeneous initial condition in the rest of the computational
domain, results in a significant difference between the two approaches in the
quantifiers that involve the gradient of $u_{P}(\bm{x},t)$. This suggests to
introduce an additional condition of continuity on the cell boundary. For zero
initial condition this can never be achieved by the strictly positive Gaussian
kernel. Therefore, we propose to alter the zero initial condition to be a
positive constant, denoted by $C$ and to investigate the selection of
$(p_{0},t_{0})$ further for such non-zero initial conditions.

In this section, we only consider the $L^{1}$-objective function from Option
$1$ and $2$ given by Equation (3.8). We take the continuity constraint for the
Gaussian-shaped extension on the boundary $\partial\Omega_{C}$, which amounts
to taking

$\displaystyle\frac{p_{0}}{4\pi Dt_{0}}\exp\\{-\frac{R^{2}}{4Dt_{0}}\\}=C$
(4.1)

as the new constraint for $(p_{0},t_{0})$ instead of Equation (3.7). By doing
this, we ensure that there exists no jump at the boundary of the cell, i.e.
$\partial\Omega_{C}$. Of course, there may be a jump in flux. Note that for
the options listed in Table 3.1 that do not include constraint, the value of
$(p_{0},t_{0})$ does not change. In Table 4.1, we present the value of
$(p_{0},t_{0})$, for Option 2 with the new constraint (4.1) with
$C=0.1,10,100$, respectively.

Table 4.1: Computed $(p_{0},t_{0})$ from the optimization problem with indicated objective function and value-pair constraint. Initial condition is $u_{0}=C$ in $\Omega\backslash\Omega_{C}$, with $C\in\\{0.1,10,100\\}$. Optimization was performed by the function optimize.minimize of Scipy package (version $1.8.0$) in Python. C |  Objective function |  Constraints |  Value of $p_{0}$ |  Value of $t_{0}$
---|---|---|---|---
$0.1$ |  $\int_{0}^{T}|\phi_{sum}(t)-\phi|dt$ |  $\displaystyle 0.1=\frac{p_{0}}{4\pi Dt_{0}}\exp\\{-\frac{R^{2}}{4Dt_{0}}\\}$ |  $4.688$ |  $0.107$
$10$ |  |  $\displaystyle 10=\frac{p_{0}}{4\pi Dt_{0}}\exp\\{-\frac{R^{2}}{4Dt_{0}}\\}$ |  $53.422$ |  $3.568$
$100$ |  |  $\displaystyle 100=\frac{p_{0}}{4\pi Dt_{0}}\exp\\{-\frac{R^{2}}{4Dt_{0}}\\}$ |  $9.194\times 10^{5}$ |  $7.310\times 10^{2}$

Figure 4.1 shows the results when $C=10$ is chosen, that is, the initial
condition is $10$ over the domain of the spatial exclusion approach. Compared
with the results in Figure 3(d) when the initial condition is zero, the
$H^{1}-$norm and the $L^{2}-$norm of the gradient of the difference between
the solutions to the two approaches have reduced. The same holds for the
quantity $c^{*}(t)$. All the subplots indicate that using of Gaussian
distribution inside the cell as initial condition reduces the difference
between the solutions to the two approaches. On the other hand, we note that
Option $1$ (red curve in Figure 4.1) performs slightly better than Option 2
with the previous constraint (black curve). This can be attributed to the fact
that the pair of $(p_{0},t_{0})$ computed without constraint in Option 1
results in the value $10.150$ approximately over the boundary of the cell.
Hence, the constraint in Equation (4.1) is almost satisfied (see the value
computed in Table 4.1), while the flux difference on the cell boundary over
time is minimized at the same time. When $C=0.1$ is chosen, similar patterns
in Figure 4.2 appear compared to that the zero initial condition is used.
However, Figure 4.2(b)-(d) show a better performance when the continuity of
the cell boundary is guaranteed; see the black curve in each subfigure. Next
to it, we also conduct the simulation with $C=100$, hence, the values of
$(p_{0},t_{0})$ computed by Option 1 and 2 will still result in a significant
difference between the two approaches, with respect to all the quantities. As
a consequence, we categorize the simulations into two subgroups: (1) Option 1
and Option 2, and (2) homogeneous initial condition and Option 2 with the new
constraint; see Figure 4.3. In all the subfigures, we observe that the results
in Subgroup (1) are at least $15$ times larger that the results in Subgroup
(2), which supports the necessity to guarantee the continuity of the initial
condition, at the cell boundary. Moreover, it can be concluded that with the
current parameter values, it is not needed to use a Gaussian-shaped extension.
The homogeneous extension of $u_{0}=C$ by the same constant value in the
cell’s interior leads to solutions that are comparable in quality as those
obtained from Gaussian-shaped extensions.

(a) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(b) $\|\nabla u_{S}(\bm{x},t)-\nabla
u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(c) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{H^{1}(\Omega\backslash\Omega_{C})}$

(d) $c^{*}(t)$

Figure 4.1: Various norms of differences between the two approaches and
$c^{*}(t)$ are shown. Here, we consider non-zero initial condition $u_{0}=C$
with $C=10$. Shown are the error quantifiers for homogeneously extended
initial condition $\bar{u}_{0}=C$ (blue) and the Gaussian-shaped extension of
the initial condition with $(p_{0},t_{0})$ computed from Option 1 (red)
without value-pair constraint, Option 2 (green) with constraint (Equation
(3.7)) in Table 3.1 and Option 2 with the new continuity constraint (Equation
(4.1)) in Table 4.1.

(a) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(b) $\|\nabla u_{S}(\bm{x},t)-\nabla
u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(c) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{H^{1}(\Omega\backslash\Omega_{C})}$

(d) $c^{*}(t)$

Figure 4.2: Various norms of differences between the two approaches and
$c^{*}(t)$ are shown. Here, we consider non-zero initial condition $u_{0}=C$
with $C=0.1$. Shown are the error quantifiers for homogeneously extended
initial condition $\bar{u}_{0}=C$ (blue) and the Gaussian-shaped extension of
the initial condition with $(p_{0},t_{0})$ computed from Option 1 (red)
without value-pair constraint, Option 2 (green) with constraint (Equation
(3.7)) in Table 3.1 and Option 2 with the new continuity constraint (Equation
(4.1)) in Table 4.1.

(a) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(b) $\|\nabla u_{S}(\bm{x},t)-\nabla
u_{P}(\bm{x},t)\|_{L^{2}(\Omega\backslash\Omega_{C})}$

(c) $\|u_{S}(\bm{x},t)-u_{P}(\bm{x},t)\|_{H^{1}(\Omega\backslash\Omega_{C})}$

(d) $c^{*}(t)$

Figure 4.3: Various norms of differences between the two approaches and
$c^{*}(t)$ are shown. Here, we consider non-zero initial condition $u_{0}=C$
with $C=100$. Shown are the error quantifiers for homogeneously extended
initial condition $\bar{u}_{0}=C$ (blue) and the Gaussian-shaped extension of
the initial condition with $(p_{0},t_{0})$ computed from Option 1 (red)
without value-pair constraint, Option 2 (green) with constraint (Equation
(3.7)) in Table 3.1 and Option 2 with the new continuity constraint (Equation
(4.1)) in Table 4.1.

For different values of $C$ in $\\{0.1,10,100\\}$ there is quite a different
total amount of mass in the environment. A relative measure of error is
therefore appropriate in order to compare the different cases, i.e. initial
conditions $u_{0}=C$ vary with $D=0.1$ fixed. We propose the ratio of the
total flux difference over the cell boundary and the total mass in
$\Omega\backslash\Omega_{C}$. By the Cauchy-Schwarz Inequality one can bound
this relative quantifier by

$\displaystyle\frac{\int_{0}^{t}\int_{\partial\Omega_{C}}|\phi(\bm{x},s)-u_{P}\nabla
D\cdot n|d\Gamma ds}{\|u_{S}(\bm{x},t)\|_{L^{1}(\Omega\backslash\Omega_{C})}}$
$\displaystyle\qquad\quad\leq\
\frac{\sqrt{|\partial\Omega_{C}|}\int_{0}^{t}\|\phi(\bm{x},s)-D\nabla
u_{P}\cdot\bm{n}ds\|_{L^{2}(\partial\Omega_{C})}ds}{\|u_{S}(\bm{x},t)\|_{L^{1}(\Omega\backslash\Omega_{C})}}\
=:\ r.e(t).$ (4.2)

The latter quantity $r.e(t)$ is more convenient because of its relation to
$c^{*}(t)$.

Note that the denominator in expression (4.2) can be computed analytically as

$\|u_{S}(\bm{x},t)\|_{L^{1}(\Omega\backslash\Omega_{C})}=\bigl{|}\Omega\backslash\Omega_{C}\bigr{|}C+2\pi
R\phi t,$

where $|\Omega\backslash\Omega_{C}|$ is the area of the domain in the spatial
exclusion approach, and $\phi$ is the constant flux density over the boundary.
Figure 4.4 shows the relative errors for various value of $C$. All the
relative errors are less than $8\%$.

Figure 4.4: The relative errors defined in Equation (4.2) are shown in the
plot, with $C=0.1,10,100$, respectively. Here, we focus on Option 2 given in
Table 3.1 and Option 2 with new constraint in Table 4.1.

##  5 Conclusions and Discussion

In this paper, we investigated conditions needed to replace the spatial
exclusion model for mass emitting objects (‘cells’) with diffusion in the
environment by a point source model for the sake of computational efficiency
or theoretical convenience in such a way that solutions to the latter still
approximate well the corresponding solutions to the former in $L^{2}$\- and
$H^{1}$-norms. Such an approximation is convenient when many cells are
considered in the computational domain and is particularly interesting when
these cells will be moving. Here, we considered the simplest case of a
circular cell, fixed at its initial location, that keeps releasing the
compounds at a constant rate. This compound is not taken up by the cell, nor
is there a reaction in the environment. It stays in the environment,
diffusing, forever after release. This simplified situation has been studied
first to identify sources of error in the approximation that stem from the
replacement of a spatial (circular) objects by a point source only.

We identified two main sources of error:

  1. (1)

A systematic time delay between the solution of the spatial exclusion model
and of the point source model, where the latter is lagging behind the former.
It originates in inappropriateness of the chosen extension of the initial
condition of the spatial exclusion model on the environment to the whole
domain, on which the point source model is defined. This time delay occurs
already for a single cell in the domain. It decreases with increasing
diffusivity. The delay can be turned into a transient effect e.g. by choosing
an extension by means of a Gaussian-shaped function on the part of the domain
that corresponds to the cell’s interior.

  2. (1)

Absence of reflection in the point source model. In this model, diffusing
compounds will move simply through the parts of the domain that constitute the
cells’ interior in the spatial exclusion model. In that model, these would
have reflected on the cell boundary and hence stayed in the environment. On
the long run, after many reflections, the two solutions should become
different.

First, we provided an analytical condition such that the solutions to the two
models are consistent, i.e. that they are equal on their shared domain of
definition, the cells’ environment. Although this condition can never be
satisfied exactly, it provides clues how to minimize the deviation between the
two approaches: minimizing the deviation between the prescribed and generated
flux density over the cell boundary. Hence, the choice for the quantity
$c^{*}(t)$, introduced in Evers et al. [3] as a measure for the quality of
approximation, is motivated. It is compared to the use of $L^{2}$\- and
$H^{1}$-norms over the environmental domain, which are more expensive to
compute.

There are infinitely many ways of extending the initial condition on the
environment to the whole domain. We chose the Gaussian-shaped extension, since
it already ‘fits’ the diffusion process. Moreover, the preliminary results
showed that all norms of differences between the two approaches converged
towards a stable value. Thus, we were inspired to extend by means of the
fundamental solution of the diffusion equation, with a chosen amplitude and
variance as parameters. Numerical results show the significant improvements in
all the quantities, such as the local norm differences and $c^{*}(t)$ and the
disappearance of the time delay after a transient.

This raised the question how to select the value of the two parameters. We
proposed and compared multiple options, which all boil down to solving single-
objective optimization problems. There seems to be no “best” option. Moreover,
a good choice may depend on which quantity is preferred for error
quantification. When a non-zero initial condition is utilized, approximation
quality improves when the extended initial condition has continuity of the
flux over the boundary of the cells. This suggests to consider a multi-
objective optimization problem in future work, where deviation from this
continuity of flux in the initial condition is balanced with a measure of
error of the flux over this boundary over time. This is a topic for future
research.

The examination of the case of multiple circular cells of identical radius in
the domain indicated the effects of absence of reflection on the quality of
approximation. This effect is also present in the situation of a single cell,
but then it is less apparent. When there are multiple cells in the domain and
their mutual distance is small compared to their radius, the solutions tend to
deviate more. In this paper, we focused on removing the systematic time delay
by proper extension of the initial condition. Further investigation is
required what conditions can be imposed such that this source of error remains
within acceptable tolerance, e.g. on this minimal distance between objects in
comparison to their diameter and the diffusion constant, and the length of the
simulation interval.

For forthcoming work, there are various possibilities of interest. In addition
to investigating the interactions between multiple cells and related
conditions mentioned above, one can think of considering the prescription of
inhomogeneous flux density over the cell boundary in space and possibly also
over time. Moreover, cells may also take up compound from the environment. In
the current paper, only circular cells have been considered, of equal radius,
as an ideal study case. However, in applications, the cells will have various
shapes. This can be resolved in a two-step approach. A spatial exclusion model
with flux over such non-circular shaped boundary may first be approximated by
such a model with inhomogeneous flux density over a circular-shaped boundary.
The latter may then be approximated by a suitable point source model. This
last type of approximation will be considered in a forthcoming paper. In the
first type of approximation, we anticipate error caused by ‘curvature
effects’. The curvature of the cell boundary influences the overall diffusion
of the chemicals in case of reflection boundary conditions. In particular,
sharp corners tend to confine diffusing matter and therefore may need
particular attention, in a sense similar to ‘absence of reflection’ as source
of error mentioned above.

In summary, this paper provided practical conditions to obtain an acceptable
approximation of the solution of a spatial exclusion model by means of a point
source model. We suggested multiple options to implement conditions
numerically. We discussed how to select the two shape-parameters of a
Gaussian-shaped extension of the initial condition of the spatial exclusion
model in a practical setting. We provide insights into deviations that will be
introduced by transferring a spatial exclusion model with diffusion to a point
source model and how two maximize the consistency between the two solutions in
their common spatial domain, over the time interval of simulation.

## References

  * Alnæs et al. [2015] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells.  The fenics project version 1.5.  _¡p¿Archive of Numerical Software_ , Vol 3:¡strong¿Starting Point and Frequency: ¡/strong¿Year: 2013¡/p¿, 2015.  doi: 10.11588/ANS.2015.100.20553.  URL http://journals.ub.uni-heidelberg.de/index.php/ans/article/view/20553.
  * Evans [2010] L. C. Evans.  _Partial differential equations_ , volume 19.  American Mathematical Soc., 2010.
  * Evers et al. [2015] J. H. Evers, S. C. Hille, and A. Muntean.  Modelling with measures: Approximation of a mass-emitting object by a point source.  _Mathematical Biosciences and Engineering_ , 12(2):357–373, 2015.  doi: 10.3934/mbe.2015.12.357.  URL https://doi.org/10.3934/mbe.2015.12.357.
  * Jackson [1999] J. D. Jackson.  Classical electrodynamics, 1999.
  * Jiao et al. [2019] Y. Jiao, Y. He, and F. Jiao.  Two-dimensional simulation of motion of red blood cells with deterministic lateral displacement devices.  _Micromachines_ , 10(6):393, jun 2019.  doi: 10.3390/mi10060393.  URL https://doi.org/10.3390%2Fmi10060393.
  * Painter and Sherratt [2003] K. J. Painter and J. A. Sherratt.  Modelling the movement of interacting cell populations.  _Journal of Theoretical Biology_ , 225(3):327–339, dec 2003.  doi: 10.1016/s0022-5193(03)00258-3.  URL https://doi.org/10.1016%2Fs0022-5193%2803%2900258-3.
  * 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.
  * Peng and Vermolen [2022a] Q. Peng and F. Vermolen.  Numerical methods to compute stresses and displacements from cellular forces: Application to the contraction of tissue.  _Journal of Computational and Applied Mathematics_ , 404:113892, Apr. 2022a.  doi: 10.1016/j.cam.2021.113892.  URL https://doi.org/10.1016/j.cam.2021.113892.
  * Peng and Vermolen [2022b] Q. Peng and F. Vermolen.  Point forces in elasticity equation and their alternatives in multi dimensions.  _Mathematics and Computers in Simulation_ , 199:182–201, Sept. 2022b.  doi: 10.1016/j.matcom.2022.03.021.  URL https://doi.org/10.1016/j.matcom.2022.03.021.
  * Peng et al. [2021] Q. Peng, F. J. Vermolen, and D. Weihs.  A formalism for modelling traction forces and cell shape evolution during cell migration in various biomedical processes.  _Biomechanics and Modeling in Mechanobiology_ , Apr. 2021.  doi: 10.1007/s10237-021-01456-2.  URL https://doi.org/10.1007/s10237-021-01456-2.
  * Schwartz [1951] L. Schwartz.  _Théorie des distributions_ , volume 1–2.  Hermann, 1951.
  * Secomb et al. [2007] T. W. Secomb, B. Styp-Rekowska, and A. R. Pries.  Two-dimensional simulation of red blood cell deformation and lateral migration in microvessels.  _Annals of Biomedical Engineering_ , 35(5):755–765, mar 2007.  doi: 10.1007/s10439-007-9275-0.  URL https://doi.org/10.1007%2Fs10439-007-9275-0.
  * Stock and Baker [2009] J. Stock and M. Baker.  Chemotaxis.  In _Encyclopedia of Microbiology_ , pages 71–78. Elsevier, 2009.  doi: 10.1016/b978-012373944-5.00068-7.  URL https://doi.org/10.1016%2Fb978-012373944-5.00068-7.
  * van Kan et al. [2005] J. van Kan, A. Segal, and F. J. Vermolen.  _Numerical methods in scientific computing_.  VSSD, 2005.
  * Wadhams and Armitage [2004] G. H. Wadhams and J. P. Armitage.  Making sense of it all: bacterial chemotaxis.  _Nature Reviews Molecular Cell Biology_ , 5(12):1024–1037, dec 2004.  doi: 10.1038/nrm1524.  URL https://doi.org/10.1038%2Fnrm1524.
  * Westphal [1968] W. H. Westphal.  Mechanics of point masses and rigid bodies.  In _A Short Textbook of Physics_ , pages 6–57. Springer Berlin Heidelberg, 1968.  doi: 10.1007/978-3-642-85476-7˙2.  URL https://doi.org/10.1007%2F978-3-642-85476-7_2.

## Appendices

##  Appendix A Solution Consistency between Two Approaches

It is out of interest to compare the two approaches, regarding the consistency
of the solutions. We recall Proposition 1.1, which shows the essential and
necessary condition to obtain consistent solutions: See 1.1

###### Proof.

Equation (1.4) is obtained as part of the proof of sufficiency of the
condition $D\nabla u_{P}(\bm{x},t)\cdot\bm{n}=\phi(\bm{x},t)$ (a.e.). To that
end, we work on the weak forms (see Section 2) of the two boundary value
problems. Note that $u_{S}(\bm{x},t)$ and $u_{P}(\bm{x},t)$ are defined on
different spatial domains. We subtract the two equations in their weak forms,
which yields

$\displaystyle\int_{\Omega\backslash\bar{\Omega}_{C}}\frac{\partial
u_{S}(\bm{x},t)}{\partial t}v_{1}(\bm{x},t)-\frac{\partial
u_{P}(\bm{x},t)}{\partial
t}v_{2}(\bm{x},t)d\Omega-\int_{\Omega_{C}}\frac{\partial
u_{P}(\bm{x},t)}{\partial
t}v_{2}(\bm{x},t)d\Omega+\int_{\partial\Omega_{C}}\phi(\bm{x},t)v_{1}(\bm{x},t)d\Gamma$
$\displaystyle+\int_{\Omega\backslash\bar{\Omega}_{C}}D\nabla
u_{S}(\bm{x},t)\nabla v_{1}(\bm{x},t)-D\nabla u_{P}(\bm{x},t)\nabla
v_{2}(\bm{x},t)d\Omega-\int_{\Omega_{C}}D\nabla u_{P}(\bm{x},t)\nabla
v_{2}(\bm{x},t)d\Omega$
$\displaystyle=-\int_{\Omega}\Phi(\bm{x},t)\delta(\bm{x}-\bm{x}_{c})v_{2}(\bm{x},t)d\Omega.$

After some simplifications (applying Gauss Theorem and boundary conditions)
and substituting $\displaystyle\frac{\partial u_{P}(\bm{x},t)}{\partial
t}=D\Delta u_{P}+\Phi(\bm{x},t)\delta(\bm{x}-\bm{x}_{c})$ in $\Omega_{C}$, we
obtain

$\displaystyle\int_{\Omega\backslash\bar{\Omega}_{C}}\frac{\partial
u_{S}(\bm{x},t)}{\partial t}v_{1}(\bm{x},t)-\frac{\partial
u_{P}(\bm{x},t)}{\partial
t}v_{2}(\bm{x},t)d\Omega+D\int_{\Omega\backslash\bar{\Omega}_{C}}\nabla
u_{S}(\bm{x},t)\nabla v_{1}(\bm{x},t)-\nabla u_{P}(\bm{x},t)\nabla
v_{2}(\bm{x},t)d\Omega$ $\displaystyle=\int_{\partial\Omega_{C}}D\nabla
u_{P}(\bm{x},t)\cdot\bm{n}v_{2}(\bm{x},t)-\phi(\bm{x},t)v_{1}(\bm{x},t)d\Gamma.$

As $v_{1}(\bm{x},t)$ and $v_{2}(\bm{x},t)$ are test functions in Hilbert
space, we select
$v_{1}(\bm{x},t)=v_{2}(\bm{x},t)=w(\bm{x},t):=u_{S}(\bm{x},t)-u_{P}(\bm{x},t)$
in $\Omega\backslash\bar{\Omega}_{C}$. Then the above equation becomes

$\displaystyle\int_{\Omega\backslash\bar{\Omega}_{C}}\frac{\partial
w(\bm{x},t)}{\partial
t}w(\bm{x},t)d\Omega+\int_{\Omega\backslash\bar{\Omega}_{C}}D\|\nabla
w(\bm{x},t)\|^{2}d\Omega=\int_{\partial\Omega_{C}}(D\nabla
u_{P}(\bm{x},t)\cdot\bm{n}-\phi(\bm{x},t))w(\bm{x},t)d\Gamma$ (A.1)
$\displaystyle\Rightarrow$ $\displaystyle\frac{\partial}{\partial
t}\int_{\Omega\backslash\bar{\Omega}_{C}}\frac{1}{2}\|w(\bm{x},t)\|^{2}d\Omega+\int_{\Omega\backslash\bar{\Omega}_{C}}D\|\nabla
w(\bm{x},t)\|^{2}d\Omega=\int_{\partial\Omega_{C}}(D\nabla
u_{P}(\bm{x},t)\cdot\bm{n}-\phi(\bm{x},t))w(\bm{x},t)d\Gamma.$

The latter equation is precisely the desired expression (1.4). Integrating
both side with respect to time $t$ from $t=0$ to any definite time $t=T$ and
using that initial conditions are equal on $\Omega\setminus\Omega_{C}$ yields

$\frac{1}{2}\|w(\bm{x},t)\|^{2}_{L^{2}(\Omega\backslash\bar{\Omega}_{C})}+D\int_{0}^{T}\|\nabla
w(\bm{x},t)\|^{2}_{L^{2}(\Omega\backslash\bar{\Omega}_{C})}dt=\int_{0}^{T}\int_{\partial\Omega_{C}}(D\nabla
u_{P}(\bm{x},t)\cdot\bm{n}-\phi(\bm{x},t))w(\bm{x},t)d\Gamma dt$ (A.2)

From Equation (A.2), if $D\nabla u_{P}(\bm{x},t)\cdot\bm{n}=\phi(\bm{x},t)$
over the boundary of the spatial exclusion $\partial\Omega_{C}$, then the
right-hand side becomes zero, subsequently, Equation (A.2) only holds when
$\|w(\bm{x},t)\|^{2}_{L^{2}(\Omega\backslash\bar{\Omega}_{C})}=0=\|\nabla
w(\bm{x},t)\|^{2}_{L^{2}(\Omega\backslash\bar{\Omega}_{C})}$, which implies

$w(\bm{x},t)=u_{S}(\bm{x},t)-u_{P}(\bm{x},t)=0,$

that is, the solutions to the two approaches are consistent in
$\Omega\backslash\bar{\Omega}_{C}$.

For the other statement, we consider the weak form of both approaches in
$(WF_{S})$ and $(WF_{P})$ and let us take a test function $v_{1}(\bm{x},t)$ on
$\Omega\setminus\Omega_{C}\times\mathbb{R}_{+}$. It can be extended to a test
function on $v_{2}$ on $\Omega\times\mathbb{R}_{+}$, which we shall denote by
the same symbol if no confusion can arise. Then we obtain

$\displaystyle\int_{\Omega\backslash\bar{\Omega}_{C}}\frac{\partial
u_{S}(\bm{x},t)}{\partial
t}v(\bm{x},t)d\Omega+\int_{\Omega\backslash\bar{\Omega}_{C}}D\nabla
u_{S}(\bm{x},t)\nabla
v(\bm{x},t)d\Omega-\int_{\partial\Omega_{C}}\phi(\bm{x},t)v(\bm{x},t)d\Gamma$
$\displaystyle=$
$\displaystyle\int_{\Omega\backslash\bar{\Omega_{C}}}\frac{\partial
u_{P}(\bm{x},t)}{\partial t}v(\bm{x},t)d\Omega+\int_{\Omega_{C}}\frac{\partial
u_{P}(\bm{x},t)}{\partial
t}v(\bm{x},t)d\Omega+\int_{\Omega\backslash\bar{\Omega_{C}}}D\nabla
u_{P}(\bm{x},t)\nabla v(\bm{x},t)d\Omega$
$\displaystyle+\int_{{\Omega_{C}}}D\nabla u_{P}(\bm{x},t)\nabla
v(\bm{x},t)d\Omega-\int_{\Omega}\Phi(\bm{x},t)\delta(\bm{x}-\bm{x}_{c})v(\bm{x},t)d\Omega.$

First, we start with proving that, if $u_{S}(\bm{x},t)=u_{P}(\bm{x},t)$, then
$\phi(\bm{x},t)-D\nabla u_{P}(\bm{x},t)\cdot\bm{n}=0,\mbox{ on
$\partial\Omega_{C}$}$, where $\bm{n}$ is the unit norm vector pointing
towards the centre of $\Omega_{C}$. Hence, given that the solutions are equal,
the above equation can be simplified and yields

$\displaystyle-\int_{\partial\Omega_{C}}\phi(\bm{x},t)v(\bm{x},t)d\Gamma+\int_{\Omega}\Phi(\bm{x},t)\delta(\bm{x}-\bm{x}_{c})v(\bm{x},t)d\Omega-\int_{\Omega_{C}}\frac{\partial
u_{P}(\bm{x},t)}{\partial t}v(\bm{x},t)d\Omega$
$\displaystyle-\int_{{\Omega_{C}}}D\nabla u_{P}(\bm{x},t)\nabla
v(\bm{x},t)d\Omega=0$

Since $\Omega_{C}$ is strictly embedded in $\Omega$, the partial differential
equation in $(BVP_{P})$ also holds for $\Omega_{C}$. Furthermore, as the cell
center $\bm{x}_{c}$ is inside $\Omega_{C}$, it can be concluded that

$\int_{\Omega}\Phi(\bm{x},t)\delta(\bm{x}-\bm{x}_{c})v(\bm{x},t)d\Omega=\int_{\Omega_{C}}\Phi(\bm{x},t)\delta(\bm{x}-\bm{x}_{c})v(\bm{x},t)d\Omega=\Phi(\bm{x}_{c},t)v(\bm{x}_{c},t).$

Then, the equation is rephrased as

$\displaystyle-\int_{\partial\Omega_{C}}\phi(\bm{x},t)v(\bm{x},t)d\Gamma+\int_{\Omega}\Phi(\bm{x},t)\delta(\bm{x}-\bm{x}_{c})v(\bm{x},t)d\Omega-\int_{\Omega_{C}}D\Delta
u_{P}(\bm{x},t)v(\bm{x},t)d\Omega$
$\displaystyle+\int_{\Omega_{C}}\Phi(\bm{x},t)\delta(\bm{x},t)v(\bm{x},t)d\Omega-\int_{{\Omega_{C}}}D\nabla
u_{P}(\bm{x},t)\nabla v(\bm{x},t)d\Omega=0$ $\displaystyle\Rightarrow$
$\displaystyle-\int_{\partial\Omega_{C}}\phi(\bm{x},t)v(\bm{x},t)d\Gamma-\int_{\Omega_{C}}D\nabla\cdot(\nabla
u_{P}(\bm{x},t)v(\bm{x},t))-D\nabla u_{P}(\bm{x},t)\nabla v(\bm{x},t)d\Omega$
$\displaystyle-\int_{{\Omega_{C}}}D\nabla u_{P}(\bm{x},t)\nabla
v(\bm{x},t)d\Omega=0$ $\displaystyle\Rightarrow$
$\displaystyle\int_{\partial\Omega_{C}}(\phi(\bm{x},t)-D\nabla
u_{P}(\bm{x},t)\cdot\bm{n})v(\bm{x},t)=0.$

The last step is done by the Gaussian Theorem [2] and $\bm{n}$ is pointing
towards the centre of $\Omega_{C}$. Note that $v(\bm{x},t)$ is a test function
in $H^{1}$ space, hence, by DuBois-Raymond lemma [14], we conclude that

$\phi(\bm{x},t)-D\nabla u_{P}(\bm{x},t)\cdot\bm{n}=0,\mbox{ on
$\partial\Omega_{C}$}.$

Hence, we proved that the solutions to both approaches are consistent in the
domain $\Omega\backslash\bar{\Omega}_{C}$ if and only if the flux over the
boundary of the hole is the same, that is,

$D\nabla u_{P}(\bm{x},t)=\phi(\bm{x},t),\mbox{ over $\partial\Omega_{C}$ and
$t\geqslant 0$.}$

∎

##  Appendix B Properties of $\phi_{sum}$

Recall that in the setting of a single Dirac source at $\bm{x}_{c}$ with mass
efflux rate $\Phi(t)$ in the infinitely extended space $\mathbb{R}^{2}$ and
with initial condition equal to $p_{0}P^{D}_{t_{0}}$, the flux density over
the boundary of the hole, $\partial\Omega_{C}$, which is a circle of radius
$R$ around $\bm{x}_{c}$, is given by

$\phi_{sum}(t)=\frac{p_{0}R}{8\pi
D(t+t_{0})^{2}}\exp\left\\{-\frac{R^{2}}{4D(t+t_{0})}\right\\}+\frac{\Phi}{2\pi
R}\exp\left\\{-\frac{R^{2}}{4Dt}\right\\},$ (B.1)

when $\Phi$ is assumed constant in time. In this section we shall summarize
various properties of $\phi_{sum}(t)$. The graphs of $\phi_{sum}(R,t)$ and
that of $\phi_{1}(R,t)$, $\phi_{2}(R,t)$, i.e. the first and second term in
(B.1), respectively, are shown in Figure B.1 for $\phi(\bm{x},t)=1$,
$\Phi=\pi$, $t_{0}=4.0$ and $p_{0}$ determined by Equation (3.7). That is,
$\phi_{sum}(0)=\phi=1$.

Figure B.1: The graphs of $\phi_{1}$, $\phi_{2}$ and their summation,
$\phi_{sum}$. Here, we use $t_{0}=4.0$, $\phi=1$ such that $\Phi=\pi$ and
$p_{0}$ is computed by Equation (3.7). Hence, $\phi_{sum}(0)\approx\phi$.

It is immediately clear that

$\phi_{sum}(0)=\frac{p_{0}R}{8\pi
Dt_{0}},\qquad\lim_{t\to\infty}\phi_{sum}(t)=\frac{\Phi}{2\pi R}.$ (B.2)

Therefore, if one wants to abstain from any systematic deviation of flux –
hence between the solutions for spatial exclusion and point source approach in
view of Proposition 1.1 – then it is necessary to require

$\frac{\Phi}{2\pi R}=\phi,$ (B.3)

which condition we shall assume to be satisfied from now on.

Define

$\alpha:=\frac{R^{2}}{4D},\qquad\beta:=\frac{p_{0}}{\pi
R\phi},\qquad\gamma:=\beta\phi,$

such that

$\phi_{sum}(t)=\frac{p_{0}}{2\pi
R}\alpha\,\frac{1}{(t+t_{0})^{2}}\,e^{-\alpha/(t+t_{0})}\
+\phi\,e^{-\alpha/t}\ =\
\mbox{$\frac{1}{2}$}\alpha\gamma\,\frac{1}{(t+t_{0})^{2}}\,e^{-\alpha/(t+t_{0})}\
+\phi\,e^{-\alpha/t}.$

A straightforward computation yields

$\phi_{sum}^{\prime}(t)\ =\
-\alpha\gamma\,\frac{t+t_{0}-\frac{1}{2}\alpha}{(t+t_{0})^{4}}\,e^{-\alpha/(t+t_{0})}\
+\ \alpha\phi\,\frac{1}{t^{2}}\,e^{-\alpha/t}.$ (B.4)

Define further

$g(t):=\frac{1}{t^{2}}e^{-\alpha/t}\qquad\mbox{and}\qquad
h(t):=\frac{t-\frac{1}{2}\alpha}{t^{4}}e^{-\alpha/t}.$

Then,

$\phi^{\prime}_{sum}(t^{*})=0\quad\mbox{{\it if and only} if}\quad
g(t^{*})=\beta\,h(t^{*}+t_{0}).$ (B.5)

One readily computes that

$g^{\prime}(t)=-\frac{1}{t^{4}}e^{-\alpha/t}(2t-\alpha)\quad\mbox{and}\quad
h^{\prime}(t)=-\frac{3}{t^{6}}e^{-\alpha/t}\bigl{(}t^{2}-\alpha
t+\mbox{$\frac{1}{6}$}\alpha^{2}\bigr{)}.$

Then, $g$ is a strictly positive function on $(0,\infty)$, with $g(t)\to 0$ as
$t\to\infty$ and $t\downarrow 0$. $g$ has a maximum value at
$t=\frac{1}{2}\alpha$ and $g^{\prime}(0)=0$. Similarly, $h(0)=0$ and $h(t)\to
0$ as $t\to\infty$. Moreover, $h(t)<0$ for $0<t<\frac{1}{2}\alpha$ and
$h(t)>0$ for $t>\frac{1}{2}\alpha$. $h^{\prime}(0)=0$ and there exists
$0<t_{-}<t_{+}$ that solve $t^{2}-\alpha t+\frac{1}{6}\alpha^{2}=0$. One has
$t_{\pm}=\alpha(\frac{1}{2}\pm\frac{1}{3}\sqrt{3})$. Clearly,
$h^{\prime}(t)<0$ for $t>t_{+}$, so $h$ has a positive maximum value at
$t_{+}>\frac{1}{2}\alpha$ and a negative minimum value at
$t_{-}<\frac{1}{2}{\alpha}$.

We obtain:

######  Lemma B.1.

For $t_{0}>0$, $\phi_{sum}$ has either zero, one, two or three critical
points, where – generically – $\phi_{sum}$ changes sign. $\phi_{sum}(t)<\phi$
for $t$ sufficiently large. Case-by-case:

  1. (i)

If there are no critical points, then $\phi_{sum}$ has a minimum at $t=0$.
Necessarily, $\phi_{sum}(0)<0$ and $\phi^{\prime}_{sum}(0)>0$.

  2. (ii)

If there is one critical point, then $\phi_{sum}$ has a minimum at this point
$t=t_{*}>0$. Necessarily, $\phi_{sum}(t_{*})<\phi$ and
$\phi^{\prime}_{sum}(0)<0$. $\phi_{sum}$ has a maximum at the boundary $t=0$.

  3. (iii)

If there are two critical point, then $\phi_{sum}$ has a minimum at some
$t_{*}>0$ with $\phi(t_{*})<\phi$ and a maximum at $0<t^{*}<t_{*}$.
Necessarily, $\phi^{\prime}_{sum}(0)>0$ and $\phi_{sum}$ has a minimum at the
boundary point $t=0$.

  4. (iv)

If there are three critical point, then $\phi_{sum}$ has a minimum at some
$t_{*,2}>0$ with $\phi(t_{*,2})<\phi$ and at $0<t_{*,1}<t_{*,2}$. Moreover,
there is a maximum at $t^{*}$ with $t_{*,1}<t^{*}<t_{2,*}$.Necessarily,
$\phi^{\prime}_{sum}(0)<0$ and $\phi_{sum}$ has a maximum at the boundary
point $t=0$.

###### Proof.

The claim on the number of critical points follows from the qualitative
properties of $g(t)$ and $h(t)$ described above and the critical point
characterisation in Equation (B.5). The intersection of $g(t)$ and $\beta
h(t+t_{0})$ are transversal, generically. Therefore, $\phi^{\prime}_{sum}(t)$
will change sign. One has

$\phi^{\prime}_{sum}(t)\ =\ -\alpha\gamma h(t+t_{0})+\alpha\phi g(t)\ =\
\alpha\phi g(t)\left[1-\beta\,\frac{h(t+t_{0})}{g(t)}\right].$

Moreover,

$\frac{h(t+t_{0})}{g(t)}\ =\
\frac{t^{2}(t+t_{0}-\frac{1}{2}\alpha)}{(t+t_{0})^{4}}\,\exp\left(\frac{\alpha
t_{0}}{t(t+t_{0}}\right)\ \to 0\qquad\mbox{as}\ t\to\infty.$

Hence, $\phi^{\prime}_{sum}(t)>0$ for $t$ sufficiently large. Since
$\phi_{sum}(t)\to\phi$ as $t\to\infty$, one must have $\phi_{sum}(t)<\phi$ for
$t$ large. In case $({\it i})$ the boundary point must then be a minimum with
$\phi_{sum}(0)<\phi$. In the other cases there must exist a minimum at a
largest $t$-value $t_{*}$ or $t_{*,2}>0$ with $\phi_{sum}(t_{*,2})<\phi$ ∎

Note that Figure B.1 shows an example of case (iv). It can be seen by careful
inspection of the graph that indeed, $\phi_{sum}(t)<\phi=1$ for $t$
sufficiently large.

##  Appendix C Locations of multiple cells in Section 3.3

In Section 3.3, the impact of the influence from other cells was discussed.
Figure C.1 show the locations of the cells corresponding to the cases in
Figure 3.2.

(a) One cell

(b) Two cell with smaller distance

(c) Two cell with larger distance

(d) Ten cells

Figure C.1: The cell locations on the domain corresponding to Section 3.3. For
every case, we always have the cell, of which the center is $(-3.5,-4)$ and
radius $\frac{1}{2}$. For this cell the measure $c^{*}(t)$ is computed in the
comparison in Section 3.3. To emphasize the locations, we show the schematics
of the spatial exclusion approach with a coarse mesh.