2cm2cm2cm2cm

These authors contributed equally to this work.

[1]Stanley J. Miklavcic These authors contributed equally to this work.

[1]Phenomics and Bioinformatics Research Centre, UniSA STEM, University of
South Australia, Adelaide, 5095, SA, Australia

2]School of Mathematics and Statistics, University of Melbourne, Parkville,
3010, VIC, Australia

# Nanoparticle uptake by a semi-permeable, spherical cell from an external
planar diffusive field. II. Numerical study of temporal and spatial
development validated using FEM.

Sandeep Santhosh Kumar<EMAIL_ADDRESS><EMAIL_ADDRESS>*  [

###### Abstract

In this paper we present a mathematical study of particle diffusion inside and
outside a spherical biological cell that has been exposed on one side to a
propagating planar diffusive front. The media inside and outside the spherical
cell are differentiated by their respective diffusion constants. A closed
form, large-time, asymptotic solution is derived by the combined means of
Laplace transform, separation of variables and asymptotic series development.
The solution process is assisted by means of an effective far-field boundary
condition, which is instrumental in resolving the conflict of planar and
spherical geometries. The focus of the paper is on a numerical comparison to
determine the accuracy of the asymptotic solution relative to a fully
numerical solution obtained using the finite element method. The asymptotic
solution is shown to be highly effective in capturing the dynamic behaviour of
the system, both internal and external to the cell, under a range of diffusive
conditions.

###### keywords:

diffusion, biological cell, Laplace transform, asymptotic approximation,
finite element method, numerical simulation

August 28, 2024

##  1 Introduction

The phenomenon of diffusion has its classical origin in the 3D movement of
particles in liquids or gases driven by a concentration gradient [1]. More
recently, the principles underlying diffusion manifest themselves in other
applied fields, such as sociology and economics, as well as in other
circumstances in chemistry and physics. A practical example in biotechnology
where diffusion is a key driving process is in molecular communication in
which organic molecules become the carriers of information between cells [2].
Another example application may also be highlighted: with the advent of
nanotechnology, engineered nanoparticles are used as drug delivery agents [3].

Accompanying these diverse applications is an increased number of complex,
diffusive scenarios. The specific complexity of interest arises in cell-cell
and particle-cell communication: a cell being exposed to the diffusing agents
from one direction. In a recent review, Ashraf et al. [4] pointed out that
different outcomes could result from different analysis methods for the same
system. In the same vein, different conclusions are possible, about a cell’s
propensity to take up particles from its environment, depending on the nature
of exposure. Circumstantial support for this can be gleaned from experimental
studies featuring a directional dependence of exposure [5]. The more
frequently encountered situation, in nature and likely also in the laboratory,
is arguably that of unidirectional exposure, which calls into question the
relevance of many model calculations. Early studies [1, 6, 7, 8], adopted the
simplifying assumption of spherical symmetry. This simplifying assumption, if
used in modern experimental situations would inevitably lead to a
misrepresentation of a cell’s own involvement in particle uptake.

Recent modelling efforts to quantify (nano)particle uptake by biological cells
that circumvent the spherical symmetry assumption have varied from spatially
featureless, multiprocess kinetic models [9, 10], to fully numerical, 3D
simulation models [11, 12, 13]. Schäfer and collaborators [14, 15] expounded a
so-called transfer function model, describing diffusion within a spherical
biological cell subject to general time and position dependent boundary
conditions on its semipermeable surface. The transfer function model was based
on a Laplace transform-eigenfunction expansion method. Despite the generality
of the method ensued by the arbitrary temporally- and spatially-dependent
boundary condition applied to the cell’s surface, the method falls short of
the needs of the present problem. The specific case, addressed herein and in
[16], of a semi-permeable biological cell exposed to a field of diffusing
particles generated by a planar source located in the unbounded, exterior
domain, requires the concurrent solution of the exterior diffusion problem in
order to satisfy the continuity conditions on the cell surface.

In [16] a closed-form, approximate solution of this combined diffusion problem
was presented. The solution supported the hypothesis that the time dependence
as nanoparticle uptake would explicitly reflect the unidirectional (as opposed
to isotropic) nature of exposure. Moreover, it was shown that only through a
scaling or modulation of this time dependence did the properties of the cell
itself (its permeability, internal diffusion constant and its capacity to
absorb particles) come into play. A comparison with a corresponding isotropic
model, for the same cell’s affinity, illustrated the importance of a proper
consideration of the external condition in determining a cell’s response to
externally supplied substances.

Due to space restrictions, no assessment of the accuracy of the asymptotic
solution was included in [16]. Also absent was a detailed study of either the
internal or external time and space dependent diffusion process. We address
these two deficiencies in this paper. It will be seen that the asymptotic
expressions derived for the large-time particle concentration inside and
around the cell are remarkably accurate when compared with a full numerical
solution of the same system using the finite element method (FEM). The
qualitative and quantitative agreement gives confidence to the subsequent
application of the asymptotic solution to derived quantities of experimental
relevance.

This report of our comparative study is organized as follows. In the next
section we outline the physical problem under study and its representative
mathematical model. The solution of the mathematical problem is summarized in
Section 3, while full mathematical details are relegated to Appendices. The
solution quoted here are fully equivalent to, but expressed more succinctly
than, the explicit asymptotic expression given in [16]. The alternative format
here is arguably more advantageous for numerical purposes. Particulars about
the FEM calculations are provided in Section 3.3. In Section 4 we compare full
numerical (FEM) solutions of the problem with our numerically implemented
asymptotic solution, under a range of parametric conditions. Both sets of
results are shown as qualitative, 2D heat maps, while specific quantitative
accuracy measures of the approximation are given in accompanying tables. The
paper concludes in Section 5 with some summary remarks and suggestions for
future work.

Cell (Region $II$)$c_{II},\,D_{II}$$z_{0}$$a$$z$$x$$R$External Medium (Region
$I$) $c_{I},\,D_{I}$ParticleReservoir$C_{0}$$\mathbf{N_{surf}}$$\mathbf{N}$
Figure 1: The physical system being modelled. A reservoir of particles at
concentration $C_{0}$ lies in the infinite half-space $z\geq z_{0}$ and
produces a diffusive front (indicated schematically by the horizontal dashed
lines) directed toward down ($z<z_{0}$) in the direction
$\mathbf{N}_{surf}=(0,0,-1)$, and encounters a spherical cell of radius $a\ll
z_{0}$ centered at the origin of a Cartesian coordinate system, $(0,0,0)$. The
outward pointing unit vector normal to the biological cell is defined as
$\mathbf{N}=\widehat{\mathbf{r}}$.

##  2 Method

###  2.1 The Physical Problem

Figure 1 depicts the physical circumstances of the system we are modelling. We
follow the variable notation introduced in [16]. A single biological cell
(Region $II$) of spherical radius $a$ and internal diffusion constant $D_{II}$
is located in a liquid medium (Region $I$) characterized by a diffusion
constant $D_{I}$. A distance $z_{0}\gg a$ away from the biological cell centre
is a planar boundary to an infinite half-space reservoir of particle solutes,
held at a fixed number concentration of $C_{0}$. Initially, both the cell
interior and the cell exterior, below the planar boundary of the reservoir,
are free of particles. Over time the solutes diffuse into Region $I$,
initially as a plane diffusive front which approaches the cell (Region $II$),
as illustrated by the thin dotted lines in Figure 1. As the diffusive front
approaches the biological cell the concentration contours (the level sets of
particle concentration) bend toward the cell to eventually match the cell’s
curvature. The particles then diffuse into the cell and propagate internally
either more slowly or faster than they would in the cell exterior, depending
on the magnitude of the ratio $A=D_{I}/D_{II}$.

###  2.2 The Mathematical Model

For convenience, the biological cell is positioned at the origin
($\bm{r}=(x,y,z)=(0,0,0)$) of a 3D coordinate system. Adopting the usual
coordinate structure, the interface of the particle reservoir is located at
$z=z_{0}$, with the plane having the normal unit vector,
$\bm{N}_{surf}=(0,0,-1)$, directed into Region $I$. We denote the local
concentration of solute particles at time $t$ in the exterior medium,
$z<z_{0}$, by $c_{I}(\mathbf{r},t)$ , while the concentration interior to the
cell at time $t$ is $c_{II}(\mathbf{r},t)$.

Although more general cell properties were considered in [16], for the
purposes of this report we shall simply adopt the conditions of continuity of
particle concentration and normal particle flux at the infinitesimally thin
cell boundary, but otherwise assume arbitrary diffusion constants, $D_{I}$ and
$D_{II}$. The governing equations in the two regions are thus,

$\displaystyle\dfrac{\partial c_{I}}{\partial t}$ $\displaystyle=$
$\displaystyle
D_{I}\nabla^{2}c_{I}=D_{I}\left[\dfrac{1}{r^{2}}\dfrac{\partial}{\partial
r}\left(r^{2}\dfrac{\partial c_{I}}{\partial
r}\right)+\dfrac{1}{r^{2}\sin\phi}\dfrac{\partial}{\partial\phi}\left(\sin\phi\dfrac{\partial
c_{I}}{\partial\phi}\right)+\dfrac{1}{r^{2}\sin^{2}\phi}\dfrac{\partial^{2}c_{I}}{\partial\theta^{2}}\right],\vspace{0.85cm}$
$\displaystyle\hskip
184.9429ptr>a,\quad\phi\in[0,\textrm{p}],\quad\theta\in[0,2\textrm{p}],\quad
t>0,$ $\displaystyle\dfrac{\partial c_{{II}}}{\partial t}$ $\displaystyle=$
$\displaystyle
D_{II}\nabla^{2}c_{II}=D_{II}\left[\dfrac{1}{r^{2}}\dfrac{\partial}{\partial
r}\left(r^{2}\dfrac{\partial c_{II}}{\partial
r}\right)+\dfrac{1}{r^{2}\sin\phi}\dfrac{\partial}{\partial\phi}\left(\sin\phi\dfrac{\partial
c_{II}}{\partial\phi}\right)+\dfrac{1}{r^{2}\sin^{2}\phi}\dfrac{\partial^{2}c_{II}}{\partial\theta^{2}}\right],\vspace{0.85cm}$
$\displaystyle\hskip
184.9429ptr<a,\quad\phi\in[0,\textrm{p}],\quad\theta\in[0,2\textrm{p}],\quad
t>0.$

These are complemented by initial and boundary conditions to ensure a unique
solution. The initial condition specifies that both regions are devoid of
particles at time $t=0$:

$\left\\{\begin{array}[]{l}c_{I}(\bm{r},t=0)=0,\quad\quad\text{ \ in Region
}I,\\\ c_{II}(\bm{r},t=0)=0,\quad\quad\text{ \ in Region
}II.\end{array}\right.$ (3)

One boundary condition represents the continuity of the normal components of
the solute surface fluxes across the cell boundary, $r=|\bm{r}|=a$:

$\left[\bm{N}\cdot\left(-D_{II}\bm{\nabla}c_{II}\right)\right](r=a,\phi,t)=\left[\bm{N}\cdot\left(-D_{I}\bm{\nabla}c_{I}\right)\right](r=a,\phi,t),\quad\phi\in[0,\textrm{p}],\leavevmode\nobreak\
t>0,$ (4)

where $\bm{N}=\widehat{\bm{r}}$ is the outward pointing unit normal to the
biological cell boundary $r=a$ (Figure 1). The second boundary condition
describes the condition of continuity of the concentration across the cell
surface.:

$c_{II}(r=a,\phi,t)=c_{I}(r=a,\phi,t),\quad\quad\phi\in[0,\textrm{p}],\leavevmode\nobreak\
t>0.$ (5)

The final boundary conditions reflect, on the one hand, the continual
provision of solutes at the interface with the particle reservoir, at
$z=z_{0}$ and on the other hand, the fact that far from both the cell and the
plane source, the concentration of solutes approaches zero:

$\left\\{\begin{array}[]{l}c_{I}(\bm{r},t)=C_{0},\hskip
128.0374ptz=z_{0},\quad\quad x,y\in\mathbb{R}^{2},\quad t>0,\\\
c_{I}(\bm{r},t)\longrightarrow 0,\hskip 125.19194ptz\rightarrow-\infty,\quad
x,y\in\mathbb{R}^{2},\quad t>0.\end{array}\right.$ (6)

###  2.3 Non-dimensional System of Equations

To simplify both the mathematical analysis and the numerical study to follow
we introduce the following dimensionless variables and parameters,

$\widetilde{\mathbf{r}}=\dfrac{\mathbf{r}}{a};\quad\rho=\dfrac{r}{a};\quad\tau=\dfrac{D_{II}t}{a^{2}};\quad
C_{1}=\dfrac{c_{I}}{C_{0}};\quad C_{2}=\dfrac{c_{II}}{C_{0}};\quad
b=\dfrac{z_{0}}{a};\quad A=\dfrac{D_{I}}{D_{II}}.$ (7)

Noting also that axisymmetry eliminates any dependence on the azimuthal angle,
$\theta$, the governing equations reduce to

$\displaystyle\dfrac{\partial C_{1}}{\partial\tau}$ $\displaystyle=$
$\displaystyle\left[\dfrac{1}{\rho^{2}}\dfrac{\partial}{\partial\rho}\left(\rho^{2}\dfrac{\partial
C_{1}}{\partial\rho}\right)+\dfrac{1}{\rho^{2}\sin\phi}\dfrac{\partial}{\partial\phi}\left(\sin\phi\dfrac{\partial
C_{1}}{\partial\phi}\right)\right],\quad\quad\rho>1,\vspace{0.85cm}$ (8)
$\displaystyle\dfrac{\partial C_{2}}{\partial\tau}$ $\displaystyle=$
$\displaystyle
A\left[\dfrac{1}{\rho^{2}}\dfrac{\partial}{\partial\rho}\left(\rho^{2}\dfrac{\partial
C_{2}}{\partial\rho}\right)+\dfrac{1}{\rho^{2}\sin\phi}\dfrac{\partial}{\partial\phi}\left(\sin\phi\dfrac{\partial
C_{2}}{\partial\phi}\right)\right],\quad 0<\rho<1,$ (9)

for $\phi\in[0,\textrm{p}],\tau>0$. In addition, we state the non-dimensional
versions of the continuity conditions at the cell surface,

$\left\\{\begin{array}[]{l}\dfrac{\partial
C_{2}}{\partial\rho}(\rho=1,\phi,\tau)=A\dfrac{\partial
C_{1}}{\partial\rho}(\rho=1,\phi,\tau),\\\
C_{2}(\rho=1,\phi,\tau)=C_{1}(\rho=1,\phi,\tau),\end{array}\right.\quad\quad\phi\in[0,\textrm{p}],\leavevmode\nobreak\
\tau>0.$ (10)

Finally, we have the far field conditions which, pointedly, are expressed in
Cartesian geometry,

$\left\\{\begin{array}[]{l}C_{1}(\widetilde{\bm{r}},\tau)=1,\hskip
133.72786pt\widetilde{z}=b,\\\ C_{1}(\widetilde{\bm{r}},\tau)\longrightarrow
0,\hskip
125.19194pt\widetilde{z}\rightarrow-\infty,\end{array}\right.\quad\quad(\widetilde{x},\widetilde{y})\in\mathbb{R}^{2},\quad\tau>0,$
(11)

where we have used the important dimensionless parameters, $A$ and $b$,
introduced above.

###  2.4 Re-specification of the Far-field Boundary Condition

In reducing the problem to one displaying axisymmetry (with no dependence on
azimuthal angle), the first of Eqs (11) nevertheless introduces a conflict of
spherical and planar geometries. We may resolve this conflict by noting that
_to a very good approximation_ the biological cell will not have any influence
on the concentration at points far from the cell. Hence, if the radius of the
biological cell were very much smaller than the distance between the planar
source and the cell centre, $a\ll z_{0}$, the time-dependent concentration of
diffusing particles on the surface of a virtual (_i.e._ , fictitious) sphere
of radius $R$ satisfying $z_{0}\geq R\gg a$, will be that determined by
unidirectional diffusion alone. This is likely a better approximation when
diffusion in the cell interior ($D_{II}$) is slower than diffusion in the
external medium ($D_{I}$). In short, this suggests an alternative to boundary
conditions Eqs (11).

It is easily shown that, for particle propagation from the planar source in
the direction of negative $z$ ($z\leq z_{0}$), the solution of the
dimensionless, unidirectional diffusion problem (with respect to independent
variables $t$ and $z$),

$\left\\{\begin{array}[]{ll}\dfrac{\partial
C_{1,\infty}}{\partial\tau}=A\dfrac{\partial^{2}C_{1,\infty}}{\partial\widetilde{z}^{2}},&\text{for}\quad\widetilde{z}<b,\quad\tau>0,\\\
C_{1,\infty}(\widetilde{z},\tau)=1,&\text{for}\quad\widetilde{z}=b,\quad\tau\geq
0,\\\ C_{1,\infty}(\widetilde{z},\tau)\longrightarrow 0,&\text{as
}\quad\widetilde{z}\longrightarrow-\infty,\quad\tau\geq 0,\end{array}\right.$
(12)

is [1],

$C_{1,\infty}(\widetilde{z},\tau)=\left(1-\text{erf}\left(\dfrac{(b-\widetilde{z})}{\sqrt{4A\tau}}\right)\right),\quad\quad\quad-\infty<\widetilde{z}\leq
b,\quad\tau\geq 0,$ (13)

where
$\text{erf}(y)=\dfrac{2}{\sqrt{\textrm{p}}}\int_{0}^{y}\mathrm{e}^{-x^{2}}\mathrm{d}x$
is the error function [17]. Consequently, the pair of boundary conditions in
Eq. (11) may be replaced with a single, inhomogeneous boundary condition at
the virtual sphere $\rho=X=R/a\gg 1$:

$C_{1}(\rho=X,\phi,\tau)=C_{1,\infty}(X\cos\phi,\tau),\quad\quad\quad\phi\in[0,\textrm{p}],\quad\tau>0.$
(14)

##  3 Asymptotic and Numerical Solution of the Diffusion Problem

###  3.1 Asymptotic approximation

The re-defined temporal-spatial diffusion problem is amenable to a Laplace
transform / separation-of-variables solution approach. With a Laplace
transform with respect to the time variable defined as [18],

$\overline{C}_{\varsigma}(\rho,\phi:p)=\mathfrak{L}\left\\{C_{\varsigma}\right\\}=\int_{0}^{\infty}\mathrm{e}^{-p\tau}C_{\varsigma}(\rho,\phi,\tau)\mathrm{d}\tau,\quad\quad\varsigma=1,2,$
(15)

the governing equations and boundary conditions become, for
$\phi\in[0,\textrm{p}]$,

$\left\\{\begin{array}[]{ll}p\overline{C}_{1}=\left[\dfrac{1}{\rho^{2}}\dfrac{\partial}{\partial\rho}\left(\rho^{2}\dfrac{\partial\overline{C}_{1}}{\partial\rho}\right)+\dfrac{1}{\rho^{2}\sin\phi}\dfrac{\partial}{\partial\phi}\left(\sin\phi\dfrac{\partial\overline{C}_{1}}{\partial\phi}\right)\right],&1<\rho<X,\vspace{0.65cm}\\\
\dfrac{p}{A}\overline{C}_{2}=\left[\dfrac{1}{\rho^{2}}\dfrac{\partial}{\partial\rho}\left(\rho^{2}\dfrac{\partial\overline{C}_{2}}{\partial\rho}\right)+\dfrac{1}{\rho^{2}\sin\phi}\dfrac{\partial}{\partial\phi}\left(\sin\phi\dfrac{\partial\overline{C}_{2}}{\partial\phi}\right)\right],&0<\rho<1,\vspace{0.5cm}\\\
\overline{C}_{1}(X,\phi,p)=\overline{C}_{1,\infty}(X\cos\phi,p),\vspace{0.2cm}&\\\
\dfrac{\partial\overline{C}_{2}}{\partial\rho}(1,\phi,p)=A\dfrac{\partial\overline{C}_{1}}{\partial\rho}(1,\phi,p),&\vspace{0.2cm}\\\
\overline{C}_{2}(1,\phi,p)=\overline{C}_{1}(1,\phi,p).&\end{array}\right.$
(16)

In Eq. (16), $p$ is the Laplace parameter corresponding to dimensionless time
$\tau$, and we have utilized the initial conditions, Eq. (3) in the transform
of the time derivative.

It may readily be confirmed that the separation-of-variables solution of Eq.
(16) is given by

$\left\\{\begin{array}[]{ll}\overline{C}_{1}(\rho,\phi,p)=\sum_{n=0}^{\infty}\left[A_{n}i_{n}\left(\rho\sqrt{p/A}\right)+B_{n}k_{n}\left(\rho\sqrt{p/A}\right)\right]P_{n}(\cos\phi),&1\leq\rho\leq
X,\\\
\overline{C}_{2}(\rho,\phi,p)=\sum_{n=0}^{\infty}D_{n}i_{n}\left(\rho\sqrt{p}\right)P_{n}(\cos\phi),&0\leq\rho\leq
1,\\\ \end{array}\right.\\\ $ (17)

for $\phi\in[0,\textrm{p}]$.

In Eq. (LABEL:TrnsfSln) $i_{n}$ and $k_{n}$ are the modified spherical Bessel
functions of the first and second kind, of order $n$ [17, 19],

$i_{n}(z)=\sqrt{\dfrac{\textrm{p}}{2z}}I_{n+1/2}(z)\quad;\quad\quad\quad
k_{n}(z)=\sqrt{\dfrac{2}{\textrm{p}z}}K_{n+1/2}(z),$ (18)

with $I_{n+1/2}$ and $K_{n+1/2}$ being the fractional order, modified Bessel
functions of first and second kind. Furthermore, to obtain a finite solution
defined on the closed interval $[0,\textrm{p}]$ we retain only
$P_{n}(\cos\phi)$, the Legendre polynomials of the first kind [17].

The expansion coefficients, $A_{n}$, $B_{n}$, and $D_{n}$, are chosen to
satisfy the Laplace-transformed boundary conditions in Eq (16). Upon
substituting the separation-of-variables solution, Eq. (LABEL:TrnsfSln), into
Eq. (16) and invoking orthogonality of Legendre functions we obtain a $3\times
3$ linear algebraic system of equations which may be readily solved. The
details of this process, while not complex, are lengthy and so are relegated
to Appendices. Some explicit results are also given in [16].

In order to return to the temporal domain, the solution in the Laplace
transformed domain must be inverted. This is achieved through evaluating
complex contour integrals of the Bromwich-type [18],

$C_{\varsigma}(\rho,\phi,\tau)=\mathfrak{L}^{-1}\left\\{\overline{C}_{\varsigma}\right\\}=\dfrac{1}{2\textrm{p}i}\int_{\sigma-i\infty}^{\sigma+i\infty}\overline{C}_{\varsigma}(\rho,\phi,p)\mathrm{e}^{p\tau}\mathrm{d}p,\quad\quad\varsigma=1,2.$
(19)

Inserting our separation-of-variables solution,
$\overline{C}_{1}(\rho,\phi,p)$ and $\overline{C}_{2}(\rho,\phi,p)$, into Eq.
(19) we obtain, at least formally, the eigenfunction expansions,

$\displaystyle C_{1}(\rho,\phi,\tau)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\mathfrak{L}^{-1}\left\\{A_{n}i_{n}\left(\rho\sqrt{p/A}\right)+B_{n}k_{n}\left(\rho\sqrt{p/A}\right)\right\\}P_{n}(\cos\phi),\quad
1\leq\rho\leq X.$ (20) $\displaystyle C_{2}(\rho,\phi,\tau)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\mathfrak{L}^{-1}\left\\{D_{n}i_{n}\left(\rho\sqrt{p}\right)\right\\}P_{n}(\cos\phi),\quad
0\leq\rho\leq 1.$ (21)

The above Laplace inversions are arguably impossible to achieve explicitly for
all times $\tau$ since the Laplace parameter not only appears in the arguments
of the explicit Bessel functions, it also appears deeply embedded in the
coefficients $A_{n}$, $B_{n}$ and $D_{n}$. Consequently, as a compromise
solution we seek only an asymptotic approximation valid for large times. The
derivation of this asymptotic approximation is outlined in the Appendices (see
also Section 3.2 of [16]). The final approximate solutions take the form of a
double series valid for $0\leq\rho\leq 1$ and a double series valid for
$1\leq\rho\leq X$. In each double series, the “outer” index denotes the order
of the Legendre polynomial in the eigenfunction series expansion and hence
represents the angular dependence, while the “inner” index denotes the order
of the asymptotic approximation in time $\tau$ (for a given $n$). The explicit
forms of these double series are exemplified by Eq. (27) in [16] for the time
and space dependent particle concentration in the biological cell interior.
For our purposes here, of numerical implementation, it suffices to simply
represent the Legendre eigenfunction series in the forms

$C_{1}=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\nu_{k}^{(n)}\,g(-X/\sqrt{A},k+n/2-1,\tau)P_{n}(\cos(\phi)),$
(22)

and

$C_{2}=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\psi_{k}^{(n)}\,g(-X/\sqrt{A},k+n/2-1,\tau)P_{n}(\cos(\phi)),$
(23)

where the time-and space-dependent factors, $g(x,m,\tau)$, and the constants,
$\nu_{k}^{(n)}$ and $\psi_{k}^{(n)}$, are derived in B.1 and B.2, and defined
by Eqs (B.78), (B.88), and (B.71), respectively. The $g(x,m,\tau)$ functions
involve confluent hypergeometric functions, which can be developed into
asymptotic series in the time variable.

###  3.2 Numerical implementation of asymptotic solution

Although they may be developed to arbitrary order, for the practical purpose
of numerical implementation we truncate the double series given in Eqs (22)
and (23). After some exploration we found that a good balance of computational
accuracy and computational efficiency could be achieved (see further comments
below) by considering terms up to and including $k=2$ and $n=5$, for the inner
and outer summations, respectively. Hence, the finite double summation
approximations, which are formally only valid for large times, that we’ve
implemented in our comparison with the fully numerical finite element
solution, are

$C_{1}\approx\sum_{n=0}^{5}\sum_{k=0}^{2}\nu_{k}^{(n)}\,g(-X/\sqrt{A},k+n/2-1,\tau)P_{n}(\cos(\phi)),\quad
1\leq\rho\leq X,\quad 0\leq\phi\leq\textrm{p},$ (24)

for outside the cell, and

$C_{2}\approx\sum_{n=0}^{5}\sum_{k=0}^{2}\psi_{k}^{(n)}\,g(-X/\sqrt{A},k+n/2-1,\tau)P_{n}(\cos(\phi)),\quad
0\leq\rho\leq 1,\quad 0\leq\phi\leq\textrm{p},$ (25)

for inside the cell.

Figure 2: An illustration of the finite element method discretization mesh,
here with $b=2$, $Q_{1}=6$ and $Q_{2}=4$. In our practical implementation we
have generally adopted the dimensions $Q_{1}=50$ and $Q_{2}=20$, and with
$b=X=10$. The biological cell periphery is indicated by the yellow half-
circle, and the cell exterior and interior are denoted F1 and F2,
respectively.

$\left.\dfrac{\partial
C_{1,2}}{\partial\tilde{x}}\right|_{\tilde{x}=0}=0$$\left.C_{1}\right|_{\tilde{z}=b}=1$$\left.\dfrac{\partial
C_{1}}{\partial\tilde{x}}\right|_{\tilde{x}=Q_{2}}=0$$\left.\dfrac{\partial
C_{1}}{\partial\tilde{z}}\right|_{\tilde{z}=(b-Q_{1})}=0$

For the explicit cases we examined it was found that terms of higher order
than $n=5$ in the Legendre function expansions contributed negligibly to the
approximation. The terms corresponding to $n=5$ for $k=0,1,2$ have a maximum
magnitude of $10^{-4}$ for those values of $\tau$ examined in the case of
$A>1$, and a maximum magnitude of $0.1$ for the corresponding case of $A<1$.
The first three terms of $\psi_{k}^{(n)}$ and $\nu_{k}^{(n)}$ each contribute
significantly to the approximations. If we omit the $\psi_{2}^{(n)}$ and
$\nu_{2}^{(n)}$ terms, the resulting approximations are less accurate for a
given $\tau$, but as $\tau$ increases the solutions nevertheless converge to
the same degree of approximation one finds with the $k=2$ terms included. The
process of explicitly computing the first three terms of $\psi_{k}^{(n)}$ and
$\nu_{k}^{(n)}$ is explained in Appendix B. Naturally, the convergence
behaviour and accuracy of the approximate solutions will depend on the
specific values of $X$, $A$ and $\rho$ that are assigned.

In the numerical comparison below we can only consider the two regions
relevant to the asymptotic approximations, which are summarized in Appendices
B.1 & B.2. That is, the numerical comparison pertains only to the region
defined by the inequality $0<\rho<X$.

In the figures appearing in the next section we highlight contour lines which
are curves joining points of equal concentration. These are the iso-contours.
With the help of these one can readily appreciate (visualize) the time-
dependent progress of the diffusive front.

###  3.3 Finite element solution

As a benchmark solution with which to compare our asymptotic solution, we used
the “Partial Differential Equation Toolbox” in Matlab to solve the system
described in Section 2. The `transient-axisymmetric` analysis module is used
to exploit the azimuthal symmetry of the system around the $z$-axis. Moreover,
it is sufficient to apply the finite element method to one half of the system
cross-section, here taken to be the $xz$-plane. The spatial variables for this
non-dimensional cross-section are defined as in Section 2.3,

$\tilde{x}=\frac{x}{a},\qquad\tilde{z}=\frac{z}{a}.$ (26)

Figure 2 shows a simple example of the elements generated by the PDE toolbox
in Matlab. The element edges labeled `E7` and `E8`, between $\tilde{z}=-1$ and
$\tilde{z}=+1$, represent the boundary of our biological cell. We maintain a
constant concentration condition, $C_{1}=1$, along the boundary edge `E2`,
which denotes the planar boundary generating the diffusive front. Along edges
`E1` and `E3` we impose the condition of zero normal flux, $\partial
C_{1}/\partial\tilde{x}=0$ and $\partial C_{1}/\partial\tilde{z}=0$,
respectively. The latter boundary conditions are required since the PDE
toolbox is incapable of treating infinite domain systems. Along the central,
axis of symmetry we also impose a zero gradient condition, $\partial
C_{1,2}/\partial\tilde{x}=0$. These boundary conditions and their locations
are shown in Figure 2. To best model the present infinite system, we capture
as large an external region as possible. We define parameters $Q_{1}$ (the
length of `E1`) and $Q_{2}$ (the length of `E3`) as being the length and
breadth, respectively, of a subset of the actual infinite spatial domain.
Relevant edges are labeled `E1` to `E8`, and the exterior and interior regions
are denoted `F1` and `F2`, respectively.

The FEM implementation involves a parameter $h_{m}$ which is the maximum
distance between nodes; the smaller is $h_{m}$, the greater is the resolution.

Figure 3: Non-dimensional particle distributions inside and outside the cell
predicted by the asymptotic approximation for $A=1$, and $b=X=2$.

##  4 Results and Discussion

###  4.1 Preliminary remarks

Before beginning a detailed numerical study, it is useful to establish the
qualitative bounds within which we anticipate the asymptotic approximation to
be valid. These qualitative confines of the model help to give a broader
perspective to the specific quantitative comparison to follow. With the
adopted model assumptions, and subsequent approximations made in the solution
process, we anticipate the following limitations on accuracy.

Consider first the influence of the virtual sphere boundary. For all times
_if_ $X$ is chosen such that $X\gtrsim O(1)$ then the asymptotic solution will
cease to have meaning since the biological cell’s presence will invariably
affect the particle concentration at the position of the virtual sphere’s
boundary at $\rho=X$. The concentration distribution around the latter
boundary will then be different at all times from those values assumed by the
1D planar solution, upon which the far field boundary condition is based.
However, at the same time it has to be acknowledged that since a FEM
simulation is necessarily also confined to a finite domain and hence is itself
burdened by effective/approximate outer boundary conditions, it may be
difficult to conclude from a comparison alone where and when, or even _if_ , a
deficiency exists on the part of the asymptotic solution. To avoid or at least
minimize any ambiguity, the majority of the results presented below were
obtained assuming a radii ratio of $X=10$. Nevertheless, for completeness, we
also consider a few cases for which the ratio $X$ is as little as $2$.

A second consideration is that of time. Since the mathematical exercise
produced a large-time, asymptotic solution we should anticipate that numerical
inaccuracies will emerge in a numerical implementation for short and
potentially medium times, revealed in a comparable FEM calculation. In other
words, the numerical implementation will in general likely be in error for
non-dimensional times $\tau\ll 1$ and potentially even for $\tau\approx O(1)$.
It is in this second consideration, other factors being equal, that the FEM
solution will be useful in revealing any errors. Since this latter prospect is
an important focus of this paper, we shall compare the two numerical works for
a range of non-dimensional times.

Given the first of the aforementioned two qualitative points we may ensure
greatest numerical accuracy in the asymptotic approximation by setting $X=b$,
which also conveniently reduces the extent of parameter space we are to
explore to that of two dimensions. For the most part we adopt this choice; the
only exception being our study of the effect of reducing the size of the
virtual sphere boundary. Our investigation comprises a semi-quantitative,
side-by-side comparison of temporal snapshots of particle distributions (in
the form of heat maps) produced by the two methods. On a more quantitative
level, for each graphical comparison we include an accompanying table of
$L_{\infty}$, $L_{1}$ and $L_{2}$ measures 111Relative as opposed to absolute
measures could also be considered. However, as particle concentrations were of
order 1 in the cases quantified (except at very short times) no greater
appreciation of accuracy is achieved by including both sets of measures.:

$\left\\{\begin{array}[]{l}L_{\infty}=\sup\left\\{\left|C_{asym}-C_{\text{\tiny
FEM}}\right|\right\\},\\\
L_{1}=\dfrac{1}{m}\sum_{m}\left|C_{asym}-C_{\text{\tiny FEM}}\right|,\\\
L_{2}=\sqrt{\dfrac{1}{m}\sum_{m}\left|C_{asym}-C_{\text{\tiny
FEM}}\right|^{2}},\end{array}\right.$ (27)

where $C_{\text{\tiny FEM}}$ and $C_{asym}$ denote the (dimensionless)
concentrations predicted by the FEM and asymptotic solution methods,
respectively. The indicated supremum and summations in Eq (27) extend either
over the biological cell interior, or the biological cell exterior (but within
the virtual sphere of radius $\rho=X$). The number of discrete points sampled
inside the biological cell was set at $m=50$, while external to the cell and
within the virtual sphere $\rho=X$, $m=4472$ points were sampled unless
otherwise indicated. Invariably, the non-dimensional spatial domain simulated
in the FEM calculations was a rectangle of size $Q_{1}\times Q_{2}=50\times
20$ containing one half of the total, mirror-symmetric region, with a non-
dimensional grid dimension (i.e., the maximum distance between FEM nodes) of
$h_{m}=0.2$.

Figure 4: Comparison of non-dimensional particle distributions inside and
outside the biological cell predicted by the asymptotic approximation (left)
and FEM simulation (right) for $A=50$, and $b=X=10$, at a non-dimensional time
of $\tau=0.1$. The FEM simulation region is a rectangle of dimensions
$Q_{1}\times Q_{2}=50\times 20$ with a grid spacing of $h_{m}=0.2$.

###  4.2 Numerical comparison

As a convenient first test case, capitalizing on the assumed continuity of
both concentration and normal flux across the membrane, we set the diffusion
ratio to unity $A=1$ so that the rates of diffusion inside and outside the
biological cell are equal. Under these conditions we would expect the cell, in
principle, to be effectively transparent to the planar front. Consequently, we
would expect the two predictions of particle concentration (FEM and
asymptotic) to be consistent with the trivial case of the biological cell
being completely absent. Thus, we would expect the particle iso-contours
(contours joining points of equal concentration) to be parallel to each other
and moreover parallel to the plane source of particles that produced them.
These iso-contours would hence travel in the direction perpendicular to the
planar source. These expectations are indeed realized as can be seen from the
contours shown in Figure 3 for the particle distributions at three non-
dimensional times. Incidentally, as a further test of the model’s resilience
we have here set the virtual sphere boundary to $X=b=2$. The FEM solution,
which is not shown, is graphically indistinguishable from the asymptotic
results shown in Figure 3.

For a diffusion constant ratio of $A=50$ (diffusion inside the biological cell
is 50 times slower than diffusion outside the cell), we show in Figures 4-5 a
sequence of snap-shots of these distributions from short to large
dimensionless times. Confirming the second of our two qualitative
expectations, we see that at $\tau=0.1$ (Figure 4) there is an obvious
discrepancy between the two results, with the asymptotic theory falsely
predicting an accumulation of particles inside the cell, despite the absence
of particles in the cell exterior. However, by $\tau=0.25$ (data not shown),
the prediction is already qualitatively in line with the FEM results, both
inside and outside the biological cell. The agreement naturally improves as
$\tau$ increases (Figure 5): from $\tau=0.75$ to $\tau=3$ in which interval
the asymptotic solution is all but identical to the FEM solution. The improved
agreement is confirmed quantitatively in Table 1 which shows a $L_{\infty}$
error of $1.9\times 10^{-4}$ inside the cell and a $L_{\infty}$ error of
$1.6\times 10^{-4}$ outside the cell at $\tau=3$. Both $L_{1}$ and $L_{2}$
error measures are comparable (to each other) and suggest, when compared with
the $L_{\infty}$ measure, that the maximum errors highlighted by $L_{\infty}$
are highly localized.

$\tau$ | 0.1 | 0.25 | 0.5 | 0.75 | 1 | 3
---|---|---|---|---|---|---
$L_{1}(\rho<1)$ | 1.79$\times 10^{-2}$ | 6.25$\times 10^{-3}$ | 2.73$\times 10^{-3}$ | 8.51$\times 10^{-5}$ | 3.18$\times 10^{-4}$ | 7.60$\times 10^{-5}$
$L_{2}(\rho<1)$ | 2.32$\times 10^{-2}$ | 8.20$\times 10^{-3}$ | 3.54$\times 10^{-3}$ | 1.24$\times 10^{-4}$ | 4.06$\times 10^{-4}$ | 8.30$\times 10^{-5}$
$L_{\infty}(\rho<1)$ | 4.50$\times 10^{-2}$ | 1.61$\times 10^{-2}$ | 6.87$\times 10^{-3}$ | 3.49$\times 10^{-4}$ | 9.30$\times 10^{-4}$ | 1.93$\times 10^{-4}$
$L_{1}(1<\rho<X)$ | 2.44$\times 10^{-2}$ | 4.14$\times 10^{-3}$ | 7.09$\times 10^{-4}$ | 2.67$\times 10^{-4}$ | 1.56$\times 10^{-4}$ | 6.66$\times 10^{-5}$
$L_{2}(1<\rho<X)$ | 5.13$\times 10^{-2}$ | 6.99$\times 10^{-3}$ | 1.18$\times 10^{-3}$ | 4.81$\times 10^{-4}$ | 2.61$\times 10^{-4}$ | 8.23$\times 10^{-5}$
$L_{\infty}(1<\rho<X)$ | 2.65$\times 10^{-1}$ | 3.32$\times 10^{-2}$ | 5.44$\times 10^{-3}$ | 2.38$\times 10^{-3}$ | 1.15$\times 10^{-3}$ | 1.63$\times 10^{-4}$
Table 1: Mean error measures ($L_{1}$, $L_{2}$, and $L_{\infty}$) as a
function of non-dimensional time $\tau$, inside and outside the cell. The
calculations were based on the case $A=50$, $b=X=10$, and for FEM grid
dimensions of $Q_{1}\times Q_{2}=50\times 20$ and a grid spacing of $h=0.2$.

The slower rate of diffusion inside the biological cell compared with
diffusion outside ($A=50\Rightarrow D_{II}=D_{I}/50$) leads to an initial
build-up in the particle concentration in the anterior half of the cell
($\phi\approx 0$) (this is certainly the case inside, but also outside,
although to a lesser extent) so those profiles or iso-contours appear to lead
those outside (i.e., are further along in $\tilde{z}$ (see, for example, the
$C=0.3$ iso-contour at $\tau=0.75$ or the $C=0.4$ iso-contour at $\tau=1$ in
Figure 5). However, in the posterior half of the cell
($\phi\approx\textrm{p}$), both inside and outside, the roles are reversed,
with the external iso-contours leading (in $\tilde{z}$) the corresponding
internal iso-contours. This is evidenced by the $C=0.2$ concentration iso-
contour at $\tau=0.75$ or the $C=0.25$ iso-contour at $\tau=1$ in Figure 5.
The external concentration gradient downstream and to the side of the cell,
possesses a significant vector component in the direction toward the axis of
symmetry, driving particles toward the cell’s shadow region. With time and
distance (downstream of the cell) the profile cross-sections tend once again
to be linear, corresponding to the asymptotic, unidirectional planar front
(see the $C=0.2$ contour at $\tau=1$, middle pair of panels in Figure 5).
These behaviorial details are captured by both the asymptotic approximation
and the FEM simulation. It is noteworthy that the tendency for the iso-
contours to regain the planar form is consistent with and validates the use of
the effective virtual sphere boundary condition at $X$, which takes on the
value of the unidirectional solution, Eq (13), for all
$\phi\in[0,\textrm{p}]$.

Figure 5: Comparison of non-dimensional particle distributions inside and
outside the cell predicted by the asymptotic approximation (left) and FEM
simulation (right) for $A=50$, $b=X=10$, and, top to bottom, at non-
dimensional times of $\tau=0.75,1,\text{ and }3$. The FEM simulation cell is a
rectangle of dimensions $Q_{1}\times Q_{2}=50\times 20$ with a grid spacing of
$h_{m}=0.2$.

Figure 6: Comparison of non-dimensional particle distributions inside and
outside the cell predicted by the asymptotic approximation (left) and FEM
simulation (right) for $A=1/50$, $b=X=10$, and, top to bottom, at non-
dimensional times of $\tau=1250,2500,\text{ and }7500$. The FEM simulation
cell is a rectangle of dimensions $Q_{1}\times Q_{2}=50\times 20$ with a grid
spacing of $h_{m}=0.2$.

For the case $A=1/50$ ($D_{II}=50D_{I}$) diffusion inside the biological cell
is faster than diffusion outside. Consequently, the relationship between iso-
contour profiles inside and outside is somewhat the reverse of the case found
with $A=50$, but is not a mirror reversal. Nevertheless, as we see from Figure
6, the more rapid particle movement inside the cell results in a depletion
zone in the cell anterior region ($\phi\approx 0$) and a build-up of particle
concentration in the posterior region ($\phi\approx\textrm{p}$), by those
particles that have been more rapidly transported through the cell. This more
rapid internal movement is responsible for an accumulation of particles in the
cell’s shadow region. This higher level of accumulation generates a
concentration gradient with a significant vector component directed away from
the symmetry axis, which of course tends to equilibrate concentrations
everywhere; the iso-contours some cell radii further downstream, again tend to
be planar. However, this occurs at a slower rate overall compared with the
$A=50$ example. In the cases shown, the diffusive propagation inside the cell
is so rapid that no contours of the representative concentration increments
are included.

In this case of internal diffusion being 50 times faster than external
diffusion, agreement between the asymptotic theory and the FEM simulation
still occurs but at significantly larger dimensionless times. The non-
dimensional times required to obtain the same level of agreement as in the
preceding case of $A=50$ are here two to four orders of magnitude greater that
unity (compare the times of the entries with roughly equivalent accuracies
($L_{1}$ and $L_{2}$ measures) in Table 1 and Table 2).

The significant difference in time orders is somewhat deceptive. In our
numerical simulation we had followed the lead of Philip [7] and Mild [8] and
defined a dimensionless time based on the diffusion constant of the cell’s
interior, $\tau=D_{II}t/a^{2}$. Consequently, the orders of magnitude
difference in non-dimensional time may be explained by noting that for the
same physical time, $t$, other things being equal, the ratio of the
dimensionless times in the two cases studied give

$\dfrac{\tau_{A>1}}{\tau_{A<1}}=\dfrac{\dfrac{0.02D_{1}t}{a^{2}}}{\dfrac{50D_{1}t}{a^{2}}}=\dfrac{1}{2500}.$

Hence, in terms of non-dimensional times, scaled as we have done, we would
expect a three order of magnitude large $\tau$ disparity for the same $t$
value. Equivalently stated, from Tables 1 and 2 we see, for example, almost
identical accuracy was achieved in the outer region at $\tau=0.25$ in the case
of $A=50$, as at $\tau=625$ in the case of $A=1/50$. These levels of accuracy
occurred approximately at

$t_{A>1}=\dfrac{\tau_{A>1}a^{2}}{D_{I}/A}\approx\dfrac{0.25a^{2}}{D_{I}/50}=12.5\dfrac{a^{2}}{D_{I}},$

in the case of $A=50$, and

$t_{A<1}=\dfrac{\tau_{A<1}a^{2}}{D_{I}/A}\approx\dfrac{625a^{2}}{50D_{I}}=12.5\dfrac{a^{2}}{D_{I}},$

in the case of $A=1/50$. That is, the same level of agreement occurring at
widely different non-dimensional times actually occur at the same real time.

$\tau$ | 250 | 625 | 1250 | 1875 | 2500 | 7500
---|---|---|---|---|---|---
$L_{1}(\rho<1)$ | 8.93$\times 10^{-6}$ | 2.51$\times 10^{-6}$ | 5.20$\times 10^{-6}$ | 3.02$\times 10^{-6}$ | 1.73$\times 10^{-5}$ | 8.64$\times 10^{-5}$
$L_{2}(\rho<1)$ | 1.01$\times 10^{-5}$ | 3.03$\times 10^{-6}$ | 5.48$\times 10^{-6}$ | 3.53$\times 10^{-6}$ | 5.75$\times 10^{-5}$ | 8.64$\times 10^{-5}$
$L_{\infty}(\rho<1)$ | 1.75$\times 10^{-6}$ | 8.80$\times 10^{-6}$ | 1.03$\times 10^{-5}$ | 1.07$\times 10^{-5}$ | 2.57$\times 10^{-5}$ | 9.38$\times 10^{-5}$
$L_{1}(1<\rho<X)$ | 2.44$\times 10^{-2}$ | 4.14$\times 10^{-3}$ | 6.68$\times 10^{-4}$ | 2.30$\times 10^{-4}$ | 1.85$\times 10^{-4}$ | 1.30$\times 10^{-4}$
$L_{2}(1<\rho<X)$ | 5.13$\times 10^{-2}$ | 6.99$\times 10^{-3}$ | 1.14$\times 10^{-3}$ | 4.34$\times 10^{-4}$ | 2.58$\times 10^{-4}$ | 1.72$\times 10^{-4}$
$L_{\infty}(1<\rho<X)$ | 2.65$\times 10^{-1}$ | 3.33$\times 10^{-2}$ | 5.44$\times 10^{-3}$ | 2.38$\times 10^{-3}$ | 1.15$\times 10^{-3}$ | 4.41$\times 10^{-4}$
Table 2: Mean error measures ($L_{1}$, $L_{2}$, and $L_{\infty}$) as a
function of non-dimensional time $\tau$, inside and outside the cell. The
calculations were based on the case $A=1/50$, $b=X=10$, and for FEM grid
dimensions of $Q_{1}\times Q_{2}=50\times 20$ and a grid spacing of $h=0.2$.

Finally, to return to the first qualitative point raised at the beginning of
this section, we investigate the effect of reducing the radius of our
effective, virtual sphere boundary, $X$.

Keeping fixed the distance to the physical source ($b=10$) the effect on the
particle distribution at $\tau=0.75$ for $A=50$, of decreasing boundary
radius, $X$ = 5, 2.5, 2 is demonstrated in Figure 7. In theory, the most
accurate example with this parameter set assumes $X=10$, this (along with the
benchmark FEM results) was shown in Figure 5. As we had chosen to maintain the
figure scale for comparison purposes, the cases for which $X<3$ highlight the
absence of perturbation solution data outside the virtual boundary. At the
boundary, $\rho=X$, and beyond, the presupposition is for the particle
distribution to be that of the unidirectional diffusion, Eq (13). A
qualitative indicator of support for the validity of this assumption is a zero
slope of the iso-contours there. From Figure 5 this qualitative support is
present for $X=10$. Now in Figure 7 it is also apparently the case for $X=5$.
However, for $X=2.5$ the iso-contours already show a nonzero slope at the
virtual sphere boundary indicating that although the concentration values at
$\rho=X$ are those of the unidirectional solution, there is a discontinuity in
the slope in the $z$-direction. Consequently, for this radius at least there
is a failure to merge smoothly with the unidirectional solution.
Interestingly, we see from Table 3, that the accuracy of the asymptotic
solution, relative to the FEM solution (which is independent of X), does not
improve significantly for virtual sphere radii larger than $X=5$, indicating
that it is not necessary to adopt the greatest possible sphere size to achieve
reasonable accuracy.

Figure 7: The effect of reducing the radius of outer, virtual sphere boundary on the non-dimensional particle distributions inside and outside the biological cell predicted by the asymptotic approximation. We show the case for $A=50$, and $b=10$, at a fixed non-dimensional time of $\tau=0.75$. The cases shown, from left to right, correspond to non-dimensional outer boundary radii of $X=5$, 2.5, and 2. The most accurate asymptotic calculation, other things being equal, as well as the FEM benchmark result, is shown in Figure 5. $X$ | 10 | 7.5 | 5 | 2.5 | 2
---|---|---|---|---|---
$m$ (exterior) | 4472 | 2495 | 1069 | 229 | 132
$L_{1}(\rho<1)$ | 8.51$\times 10^{-5}$ | 7.97$\times 10^{-5}$ | 1.16$\times 10^{-4}$ | 8.99$\times 10^{-4}$ | 1.80$\times 10^{-3}$
$L_{2}(\rho<1)$ | 1.24$\times 10^{-4}$ | 1.22$\times 10^{-4}$ | 1.59$\times 10^{-4}$ | 1.10$\times 10^{-3}$ | 2.12$\times 10^{-3}$
$L_{\infty}(\rho<1)$ | 3.49$\times 10^{-4}$ | 3.27$\times 10^{-4}$ | 4.19$\times 10^{-4}$ | 2.20$\times 10^{-3}$ | 4.20$\times 10^{-3}$
$L_{1}(1<\rho<X)$ | 2.67$\times 10^{-4}$ | 1.46$\times 10^{-4}$ | 3.60$\times 10^{-4}$ | 1.90$\times 10^{-3}$ | 3.20$\times 10^{-3}$
$L_{2}(1<\rho<X)$ | 4.81$\times 10^{-4}$ | 1.86$\times 10^{-4}$ | 4.15$\times 10^{-4}$ | 2.08$\times 10^{-3}$ | 3.53$\times 10^{-3}$
$L_{\infty}(1<\rho<X)$ | 2.38$\times 10^{-3}$ | 6.34$\times 10^{-4}$ | 8.23$\times 10^{-4}$ | 3.60$\times 10^{-3}$ | 5.70$\times 10^{-3}$
Table 3: Mean error measures ($L_{1}$, $L_{2}$, and $L_{\infty}$) inside and
outside the cell as a function of radius of outer, virtual spherical boundary
$X$. The calculations were based on the case $A=50$, $b=10$, and at a non-
dimensional time of $\tau=0.75$. The number of interior points sampled is
fixed at $m=50$.

##  5 Summary and Concluding Remarks

The asymptotic solution to the problem of particle diffusion in and around a
spherical cell exposed to a plane source of particles is compared with a fully
numerical solutions using the finite element method. We have indicated regimes
in which the asymptotic solution is expected to perform poorly or fail
altogether, as well as those regimes when it is expected to be accurate. The
arguments have been confirmed by the numerical comparison with the FEM
solutions. In fact, the asymptotic solution is shown to be surprisingly
accurate over a considerable range of system parameters, even when the
Legendre series is truncated at $n=5$ and the asymptotic series is truncated
at $k=2$.

We found that good agreement, often to four decimal places, between the
asymptotic approximation and the FEM solution occurred at considerably
different non-dimensional times for different $A$ values. However, it was also
found that these different non-dimensional times were in correspondence with
the relative values of the diffusion constants, and agreement occurred at the
same physical time.

The principal aim achieved in this work was the demonstration of good
agreement between our asymptotic solution and a conventional, finite element
solution of the diffusion equations. The agreement validates the use of the
asymptotic solution expression (explicitly presented in [16]), instead of less
convenient tables of numerical values, to determine derived quantities that
are accessible experimentally. On this last note, an explicit example
application is considered in [20].

An experimental detail, not yet explored theoretically but clearly warranted
as indicated by the experimental studies, was the influence of gravity in the
diffusion process. It was found by experiment [5] that the orientation of the
source-to-cell line relative to the direction of the gravitational force
affected the rate at which nanoparticles accumulated in cells. Thus, it
remains to include sedimentation in a future model. Such a theoretical
analysis of the combined effect of sedimentation and diffusion would ascertain
to what degree and in what form sedimentation affects the time dependent
accumulation of nanoparticles in a spherical cell. We shall pursue this
problem and will report on our findings, including a comparison with
experimental data, in a separate publication.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

SJM was responsible for the model development, mathematical methodology,
research direction, interpreting of results and drafting of the paper. SSK was
responsible for independently verifying the analysis, undertaking the
numerical study (FEM and asymptotic approximation), and preparing figures.
Both authors reviewed and approved the final version of the manuscript.

## Funding

This project is supported by the Australian Research Council (Discovery
project grant DP200103168).

## Acknowledgments

## References

  *   * Crank [1956] Crank, J.: The Mathematics of Diffusion. Oxford University Press, Oxford (1956)
  * Farsad et al. [2016] Farsad, N., Yilmaz, H.B., Eckford, A., Chae, C.-B., Guo, W.: A comprehensive survey of recent advancements in molecular communication. IEEE Communications, Surveys and Tutorials 18(3), 1887–1919 (2016)
  * Allen and Cullis [2004] Allen, T.M., Cullis, P.R.: Drug delivery systems: entering the mainstream. Science 303, 1818–1822 (2004)
  * Ashraf et al. [2020] Ashraf, S., Said, A.H., Hartmann, R., Assmann, M.-A., Feliu, N., Lenz, P., Parak, W.J.: Quantitative particle uptake by cells as analyzed by different methods. Angewandte Chemie International Edition 59, 5438–5453 (2020)
  * Cui et al. [2016] Cui, J., Faria, M., Bjornmalm, M., Ju, Y., Suma, T., Gunawan, S.T., Richardson, J.J., Heidari, H., Bals, S., Crampin, E.J., Caruso, F.: A framework to account for sedimentation and diffusion in particle-cell interactions. Langmuir 32, 12394–12402 (2016)
  * Rashevsky [1948] Rashevsky, N.: Mathematical Biophysics. University of Chicago Press, Chicago (1948)
  * Philip [1964] Philip, J.R.: Transient heat conduction between a sphere and a surrounding medium of different thermal properties. Australian Journal of Physics 17, 423–430 (1964)
  * Mild [1971] Mild, K.H.: The kinetics of diffusion between a spherical cell and a surrounding medium with different diffusion properties. Bulletin of Mathematical Biophysics 33, 17–26 (1971)
  * Sorrell et al. [2014] Sorrell, I., Shipley, R.J., Hearnden, V., Colley, H.E., Thornhill, M.H., Murdoch, C., Webb, S.D.: Combined mathematical modelling and experimentation to predict polymersome uptake by oral cancer cells. Nanomedicine: Nanotechnology, Biology and Medicine 10, 339–348 (2014)
  * West et al. [2021] West, H., Roberts, F., Sweeney, P., Walker-Samuel, S., Leedale, J., Colley, H., Murdoch, C., Shipley, R.J., Webb, S.D.: A mathematical investigation into the uptake kinetics of nanoparticles in vitro. PLoS ONE 16(7), 0254208 (2021)
  * Hinderliter et al. [2010] Hinderliter, P.M., Minard, K.R., Orr, G., Chrisler, W.B., Thrall, B.D., Pounds, J.G., Teeguarden, J.G.: Issd: A computational model of particle sedimentation, diffusion and target cell dosimetry for in vitro toxicity studies. Particle and Fibre Toxicology 7, 36–120 (2010)
  * Al-Obaidi and Florence [2015] Al-Obaidi, H., Florence, A.T.: Nanoparticle delivery and particle diffusion in confined and complex environments. Journal of Drug Delivery Science and Technology 30, 266–277 (2015)
  * Friedmann [2016] Friedmann, E.: PDE/ODE modeling and simulation to determine the role of diffusion in long-term and -range cellular signaling. BMC Biophysics 8, 10–116 (2016)
  * Schäfer and Rabenstein [2019] Schäfer, M., Rabenstein, R.: An analytical model of diffusion in a sphere with semi-permeable boundary. Proc. Workshop in Molecular Communications, 1–2 (2019)
  * Schäfer et al. [2020] Schäfer, M., Wicke, W., Haselmayr, W., Rabenstein, R., Schober, R.: Spherical diffusion model with semi-permeable boundary: a transfer function approach. IEEE International Conference on Communications ICC, 1–7 (2020)
  * Miklavcic [2024] Miklavcic, S.J.: Nanoparticle uptake by a semi-permeable spherical cell from an external, planar diffusive field. I. Mathematical model and asymptotic solution. Journal of Engineering Mathematics, (2024)
  * Abramowitz and Stegun [1964] Abramowitz, M., Stegun, I.A.: Handbook of Mathematical Functions. Dover Publications, New York (1964)
  * Morse and Feshbach [1953] Morse, P.M., Feshbach, H.: Methods of Theoretical Physics. McGraw-Hill Book Company, New York (1953)
  * Arfken and Weber [2001] Arfken, B. George, Weber, H.J.: Mathematical Methods for Physicists, 5th edn. Academic Press, San Diego (2001)
  * Miklavcic [2024] Miklavcic, S.J.: Nanoparticle uptake by a semi-permeable spherical cell from an external, planar diffusive field. III. Analysis of an experimental system. Journal of Engineering Mathematics, (2024)
  * Hahn and Özisik [2012] Hahn, D.W., Özisik, M.N.: Heat Conduction, 3rd edn. Wiley, Hoboken, New Jersey (2012)
  * Carslaw and Jaeger [1941] Carslaw, H.S., Jaeger, J.C.: Operational Methods in Applied Mathematics, 1st edn. Clarendon Press, Oxford (1941)
  * Olver et al. [2021] Olver, F.W.J., Olde Daalhuis, A.B., Lozier, D.W., Schneider, B.I., Boisvert, R.F., Clark, C.W., Miller, B.R., Saunders, B.V., Cohl, H.S., M. A. McClain, e.: NIST Digital Library of Mathematical Functions. Release 1.1.2 (2021). http://dlmf.nist.gov/
  * Weisstein [2020] Weisstein, E.W.: Confluent Hypergeometric Function of the First Kind. From MathWorld—A Wolfram Web Resource (2020). https://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html
  * Wolfram Research [2020] Wolfram Research: Inverse Laplace Transform (2020). https://reference.wolfram.com/language/ref/InverseLaplaceTransform.html

.

##  Appendix A Separation of Variables Solution in the Laplace Domain

The function $\overline{C}_{2}$,

$\overline{C}_{2}=\sum_{n=0}^{\infty}(D_{n}i_{n}(\rho\sqrt{p})+E_{n}k_{n}(\rho\sqrt{p}))P_{n}(\cos(\phi)),$
(A.28)

describes the Laplace transformed concentration inside the cell. Analogously,
the function $\overline{C}_{1}$,

$\overline{C}_{1}=\sum_{n=0}^{\infty}(A_{n}i_{n}(\rho\sqrt{p/A})+B_{n}k_{n}(\rho\sqrt{p/A}))P_{n}(\cos(\phi)),$
(A.29)

represents the Laplace transformed concentration outside the cell. To fully
determine these solutions we enforce the boundary conditions in the Laplace
domain to determine the coefficients $A_{n}$, $B_{n}$, $D_{n}$ and $E_{n}$.
For example, we observe that to eliminate the singularity at the origin
$\rho=0$ we require

$E_{n}=0,\quad n=0,1,2,\ldots\,,$ (A.30)

Furthermore, to ensure continuity of concentration across the cell surface (at
$\rho=1$) we require,

$A_{n}i_{n}(\sqrt{p/A})+B_{n}k_{n}(\sqrt{p/A})=D_{n}i_{n}(\sqrt{p})\quad
n=0,1,2,\ldots\,.$ (A.31)

Similarly, to ensure continuity of the flux density across the cell surface we
require,

$\sqrt{A}\Big{(}A_{n}i_{n}^{\prime}(\sqrt{p/A})+B_{n}k_{n}^{\prime}(\sqrt{p/A})\Big{)}=D_{n}i_{n}^{\prime}(\sqrt{p})$
(A.32)

where for brevity we define

$i_{n}^{\prime}(x)=\frac{\partial(i_{n}(x))}{\partial x},\quad\quad
k_{n}^{\prime}(x)=\frac{\partial(k_{n}(x))}{\partial x},\quad
n=0,1,2,\ldots\,.$ (A.33)

Finally, we impose the far-field condition on the virtual spherical boundary,
$X$. Utilising the expression known as the plane-wave expansion which can be
found in [17](Eq. 10.1.47), yields

$A_{n}i_{n}(X\sqrt{p/A})+B_{n}k_{n}(X\sqrt{p/A})=\frac{(2n+1)}{p}e^{-b\sqrt{p/A}}i_{n}(X\sqrt{p/A}),\quad
n=0,1,\ldots\,.$ (A.34)

Equations A.31, A.32 & A.34 form a system of three linear equations in the
unknown parameters $A_{n}$, $B_{n}$ and $D_{n}$.

To solve this system we first introduce the quantities,222Many of the
following parameters and variables should additionally carry the index “$n$”.
However, in the interest of clarity we have suppressed the $n$. The reader
should be mindful of the implicit dependence.

$\displaystyle\sigma_{l}=i_{n}(x_{l}),\quad\omega_{l}$
$\displaystyle=k_{n}(x_{l}),\quad\sigma_{l}^{\prime}=\frac{\partial(i_{n}(x_{l}))}{\partial
x_{l}},\quad\omega_{l}^{\prime}=\frac{\partial(k_{n}(x_{l}))}{\partial
x_{l}},\quad l=1,2,\ldots,5$ (A.35) where
$\displaystyle\\{x_{1},x_{2},\ldots,x_{5}\\}=\\{\sqrt{p},\,\rho\sqrt{p},\,\sqrt{\frac{p}{A}},\,\rho\sqrt{\frac{p}{A}},\,X\sqrt{\frac{p}{A}}\\}.$
(A.36)

Solving the system of linear equations in terms of these new quantities, we
obtain

$\displaystyle A_{n}$
$\displaystyle=\Big{(}\frac{(2n+1)}{p}e^{-b\sqrt{\frac{p}{A}}}\Big{)}\frac{-(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime})\sigma_{5}}{\omega_{5}(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})-\sigma_{5}(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime})},$
(A.37) $\displaystyle B_{n}$
$\displaystyle=\Big{(}\frac{(2n+1)}{p}e^{-b\sqrt{\frac{p}{A}}}\Big{)}\frac{(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})\sigma_{5}}{\omega_{5}(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})-\sigma_{5}(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime})},$
(A.38) $\displaystyle D_{n}$
$\displaystyle=\Big{(}\frac{(2n+1)}{p}e^{-b\sqrt{\frac{p}{A}}}\Big{)}\frac{\sqrt{A}(\sigma_{3}^{\prime}\omega_{3}-\sigma_{3}\omega_{3}^{\prime})\sigma_{5}}{\omega_{5}(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})-\sigma_{5}(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime})}.$
(A.39)

To summarize, the transformed solutions, $\overline{C}_{1}$ and
$\overline{C}_{2}$, for the concentration outside and inside the cell,
respectively, are

$\displaystyle\overline{C}_{2}$
$\displaystyle=\sum_{n=0}^{\infty}\overline{V}_{n}P_{n}(\cos(\phi)),$ (A.40)
where $\displaystyle\overline{V}_{n}$
$\displaystyle=\Big{(}\frac{(2n+1)}{p}e^{-b\sqrt{\frac{p}{A}}}\Big{)}\frac{\sqrt{A}(\sigma_{3}^{\prime}\omega_{3}-\sigma_{3}\omega_{3}^{\prime})\sigma_{5}\sigma_{2}}{\omega_{5}(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})-\sigma_{5}(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime})}$
(A.41)

and

$\displaystyle\overline{C}_{1}$
$\displaystyle=\sum_{n=0}^{\infty}\overline{U}_{n}P_{n}(\cos(\phi)),$ (A.42)
with $\displaystyle\overline{U}_{n}$
$\displaystyle=\Big{(}\frac{(2n+1)}{p}e^{-b\sqrt{\frac{p}{A}}}\Big{)}\frac{-(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime})\sigma_{5}\sigma_{4}+(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})\sigma_{5}\omega_{4}}{\omega_{5}(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})-\sigma_{5}(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime})}.$
(A.43)

.

##  Appendix B Asymptotic Approximation

The solutions in the Laplace domain given by Eqs. A.40–A.43 are complicated
given the appearance of the Laplace parameter in the arguments of the modified
spherical Bessel functions and in the majority of factors multiplying the
Bessel functions. Consequently, achieving a closed form expression for the
inverse Laplace transform is highly unlikely. Instead, as a compromise
solution we proceed with an asymptotic analysis.

It is known that the Laplace transform $\overline{f}(p)$ of a function as
$p\to 0$ corresponds to an approximation for $f(\tau)$ valid for
$\tau\to\infty$. This principle is discussed by Hahn and Özisik [21] and is
used by Philip [7] and Mild [8]. A similar concept is formulated by Carslaw
and Jaeger [22]. We shall take advantage of this fact to obtain asymptotic
approximations to the Laplace inverse of Eqs. A.40–A.43 valid for large times.

###  B.1 Large-time, asymptotic approximation for the concentration inside the
cell

We begin by developing a series expansion of $\overline{C}_{2}$ as $p\to 0$.
We rewrite Eqs. A.40 & A.41 in terms of our new quantities, Eqs. A.35 & A.36,

$\overline{C}_{2}=\sum_{n=0}^{\infty}\overline{V}_{n}P_{n}(\cos(\phi)),\quad\overline{V}_{n}=\Big{(}\frac{(2n+1)}{p}e^{-b\sqrt{\frac{p}{A}}}\Big{)}\frac{\Omega_{1}}{\Omega_{2}}$
(B.44)

where

$\Omega_{1}=\sqrt{A}(\sigma_{3}^{\prime}\omega_{3}-\sigma_{3}\omega_{3}^{\prime})\sigma_{5}\sigma_{2},$
(B.45)
$\Omega_{2}=\omega_{5}(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})-\sigma_{5}(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime}).$
(B.46)

We first determine an expansion for $1/\Omega_{2}$. We introduce the
quantities

$\beta_{l}=i_{-n}(x_{l}),\quad\beta_{l}^{\prime}=\frac{\partial(i_{-n}(x_{l}))}{\partial
x_{l}},\quad l=1,2,\ldots,5$ (B.47)

where $i_{-n}(z)$ is the negative order modified spherical Bessel function of
the first kind and with $x_{l}$ as defined in Eq. A.36. Using Eq. C.95 (in
C.1), we have

$\omega_{l}=\frac{(-1)^{n+1}\textrm{p}}{2}(\sigma_{l}-\beta_{l})$ (B.48)

and

$\omega_{l}^{\prime}=\frac{(-1)^{n+1}\textrm{p}}{2}(\sigma_{l}^{\prime}-\beta_{l}^{\prime}).$
(B.49)

We substitute Eqs. B.48 & B.49 into Eq. B.46, to obtain

$\Omega_{2}=\frac{(-1)^{n+1}\textrm{p}}{2}(\beta_{5}(\sqrt{A}\sigma_{1}\sigma_{3}^{\prime}-\sigma_{1}^{\prime}\sigma_{3})-\sigma_{5}(\sqrt{A}\sigma_{1}\beta_{3}^{\prime}-\sigma_{1}^{\prime}\beta_{3})).$
(B.50)

We now introduce an expansion for each $\sigma_{l}$, which are the modified
spherical Bessel functions of the first kind, using Eq. C.97:

$\sigma_{l}=(\sqrt{p})^{n}\sum_{k=0}^{\infty}\sigma_{l,k}\,(\sqrt{p})^{2k}$
(B.51)

where

$\sigma_{l,k}=\frac{a(n,k)(x_{l})^{2k+n}}{(\sqrt{p})^{2k+n}},\qquad
k=0,1,\ldots\,,$ (B.52)

$a(n,k)$, defined in Eq. C.98 are coefficients that are independent of $p$. We
note that $\sigma_{l,k}$ are the coefficients of $(\sqrt{p})^{2k+n}$ in the
expansion for $i_{n}(x_{l})$ using Equation C.97. We observe that $x_{l}$
(Equation A.36) is comprised of $\sqrt{p}$ in the numerator, for each $l$. The
factor of $\sqrt{p}$ is eliminated from $\sigma_{l,k}$ by the
$(\sqrt{p})^{2k+n}$ in the denominator of Equation B.52. Hence, $\sigma_{l,k}$
for all values of $k$ and $l$ are independent of $p$. Similarly, we introduce
expansions for $\sigma_{l}^{\prime}$, $\beta_{l}$ and $\beta_{l}^{\prime}$
using Equations C.97, C.99 & C.100

$\sigma^{\prime}_{l}=(\sqrt{p})^{(n-1)}\sum_{k=0}^{\infty}\sigma^{\prime}_{l,k}\,(\sqrt{p})^{2k},\quad\beta_{l}=(\sqrt{p})^{-(n+1)}\sum_{k=0}^{\infty}\beta_{l,k}\,(\sqrt{p})^{2k},$
(B.53)
$\beta_{l}^{\prime}=(\sqrt{p})^{-(n+2)}\sum_{k=0}^{\infty}\beta_{l,k}\,(\sqrt{p})^{2k},\quad
l=1,2,\ldots,5$ (B.54)

where

$\sigma_{l,k}^{\prime}=\frac{a(n,k)(2k+n)(x_{l})^{2k+(n-1)}}{(\sqrt{p})^{2k+(n-1)}},\quad\beta_{l,k}=\frac{b(n,k)(x_{l})^{2k-(n+1)}}{(\sqrt{p})^{2k-(n+1)}},\quad$
(B.55)
$\beta_{l,k}^{\prime}=\frac{b(n,k)(2k-n-1)(x_{l})^{2k-(n+2)}}{(\sqrt{p})^{2k-(n+2)}},\quad
k=0,1,...\,.$ (B.56)

Analogous to $\sigma_{l,k}$, we observe that $\sigma_{l,k}^{\prime}$,
$\beta_{l,k}$ and $\beta_{l,k}^{\prime}$ are not dependent on $p$ for all
values of $l$ and $k$. Substituting Eqs. B.51, B.53 & B.54 into Eq. B.50, we
obtain

$\Omega_{2}=\frac{(-1)^{n+1}\textrm{p}\,p^{n/2-1}}{2}(\Phi_{1}-\Phi_{2})$
(B.57)

where

$\displaystyle\Phi_{1}$
$\displaystyle=\Big{(}\sum_{k=0}^{\infty}\beta_{5,k}\,p^{k}\Big{)}\Bigg{(}\sqrt{A}\Big{(}\sum_{k=0}^{\infty}\sigma_{1,k}\,p^{k}\Big{)}\Big{(}\sum_{k=0}^{\infty}\sigma_{3,k}^{\prime}\,p^{k}\Big{)}-\Big{(}\sum_{k=0}^{\infty}\sigma_{1,k}^{\prime}\,p^{k}\Big{)}\Big{(}\sum_{k=0}^{\infty}\sigma_{3,k}\,p^{k}\Big{)}\Bigg{)},$
(B.58) $\displaystyle\Phi_{2}$
$\displaystyle=\Big{(}\sum_{k=0}^{\infty}\sigma_{5,k}\,p^{k}\Big{)}\Bigg{(}\sqrt{A}\Big{(}\sum_{k=0}^{\infty}\sigma_{1,k}\,p^{k}\Big{)}\Big{(}\sum_{k=0}^{\infty}\beta_{3,k}^{\prime}\,p^{k}\Big{)}-\Big{(}\sum_{k=0}^{\infty}\sigma_{1,k}^{\prime}\,p^{k}\Big{)}\Big{(}\sum_{k=0}^{\infty}\beta_{3,k}\,p^{k}\Big{)}\Bigg{)}.$
(B.59)

The expressions for $\Phi_{1}$ and $\Phi_{2}$ are composed of terms which are
products of three series in $p$. Using the triple Cauchy product formula from
Eq. C.104, we get the single series for the denominator in Eq. B.44

$\Omega_{2}=p^{n/2-1}\sum_{k=0}^{\infty}\gamma_{k}\,p^{k}$ (B.60)

where

$\gamma_{k}=\frac{(-1)^{n+1}\textrm{p}}{2}\sum_{i=0}^{k}\sum_{j=0}^{i}\beta_{5,k-i}(\sqrt{A}\sigma_{1,i-j}\,\sigma_{3,j}^{\prime}-\sigma_{1,i-j}^{\prime}\,\sigma_{3,j})-\sigma_{5,k-i}(\sqrt{A}\sigma_{1,i-j}\,\beta_{3,j}^{\prime}-\sigma_{1,i-j}^{\prime}\,\beta_{3,j}).$
(B.61)

Manipulating Eq. B.60, we obtain the expression

$\frac{1}{\Omega_{2}}=\frac{p^{(1-n/2)}}{\gamma_{0}}\Big{(}1+\sum_{k=1}^{\infty}\frac{\gamma_{k}}{\gamma_{0}}\,p^{k}\Big{)}^{-1}$
(B.62)

Using the geometric series expansion [21], the last equation may be re-written
as

$\frac{1}{\Omega_{2}}=p^{(1-n/2)}\sum_{r=0}^{\infty}\frac{(-1)^{r}}{\gamma_{0}}\Big{(}\sum_{k=1}^{\infty}\frac{\gamma_{k}}{\gamma_{0}}\,p^{k}\Big{)}^{r}=p^{(1-n/2)}\sum_{k=0}^{\infty}\mu_{k}\,p^{k},$
(B.63)

which is valid for small values of $p$ since

$\Big{|}\sum_{k=1}^{\infty}\frac{\gamma_{k}}{\gamma_{0}}p^{k}\Big{|}<1,\qquad\textrm{as}\qquad
p\to 0.$ (B.64)

We determine the first four terms of $\mu_{k}$,

$\mu_{0}=\frac{1}{\gamma_{0}},\quad\mu_{1}=-\frac{\gamma_{1}}{\gamma_{0}^{2}},\quad\mu_{2}=\frac{\gamma_{1}^{2}}{\gamma_{0}^{3}}-\frac{\gamma_{2}}{\gamma_{0}^{2}},\quad\mu_{3}=-\frac{\gamma_{1}^{3}}{\gamma_{0}^{4}}+2\frac{\gamma_{1}\gamma_{2}}{\gamma_{0}^{3}}-\frac{\gamma_{3}}{\gamma_{0}^{2}}\,.$
(B.65)

We now develop a series expansion for $\Omega_{1}$ in the numerator of Eq.
B.44. We observe that the term
$\sigma_{3}^{\prime}\omega_{3}-\sigma_{3}\omega_{3}^{\prime}$ in Eq. B.45 is
the Wronskian of $i_{n}(x_{3})$ and $k_{n}(x_{3})$ with $x_{3}=\sqrt{p/A}$.
Hence, using Eq. C.96, we find

$\sigma_{3}^{\prime}\omega_{3}-\sigma_{3}\omega_{3}^{\prime}=-\frac{\textrm{p}A}{2p}.$
(B.66)

Substitution of this into Eq. B.45, we get

$\Omega_{1}=-\frac{A^{3/2}\textrm{p}}{2p}\sigma_{5}\sigma_{2}.$ (B.67)

Using Eq. B.51, we obtain

$\Omega_{1}=-\frac{A^{3/2}\textrm{p}}{2}p^{n-1}\Big{(}\sum_{k=0}^{\infty}\sigma_{5,k}\,p^{k}\Big{)}\Big{(}\sum_{k=0}^{\infty}\sigma_{2,k}\,p^{k}\Big{)}.$
(B.68)

We substitute Eqs. B.63 & B.68 into B.44, to get our Legendre coefficient
functions,

$\overline{V}_{n}=-\frac{(2n+1)A^{3/2}\textrm{p}}{2}e^{-b\sqrt{\frac{p}{A}}}p^{n/2-1}\Big{(}\sum_{k=0}^{\infty}\mu_{k}\,p^{k}\Big{)}\Big{(}\sum_{k=0}^{\infty}\sigma_{5,k}\,p^{k}\Big{)}\Big{(}\sum_{k=0}^{\infty}\sigma_{2,k}\,p^{k}\Big{)}.$
(B.69)

Using the triple Cauchy product formula Eq. C.104, we arrive at the final
series expansion

$\overline{V}_{n}=e^{-b\sqrt{\frac{p}{A}}}p^{n/2-1}\sum_{k=0}^{\infty}\psi_{k}^{(n)}\,p^{k},$
(B.70)

where

$\psi_{k}^{(n)}=-\frac{(2n+1)A^{3/2}\textrm{p}}{2}\sum_{i=0}^{k}\sum_{j=0}^{i}\mu_{k-i}\,\sigma_{5,i-j}\,\sigma_{2,j},\quad
n=0,1,2\ldots\,.$ (B.71)

We note that the radial dependence is contained in the factor $\sigma_{2,j}$,
but more importantly, none of these factors contain $p$. We now proceed to
determine the inverse Laplace transform of $\overline{C}_{2}$, using Eqs. B.44
& B.70,

$\mathcal{L}^{-1}\\{\overline{C}_{2}\\}=\sum_{n=0}^{\infty}\mathcal{L}^{-1}\\{\overline{V}_{n}\\}P_{n}(\cos(\phi)).$
(B.72)

where

$\mathcal{L}^{-1}\\{\overline{V}_{n}\\}=\sum_{k=0}^{\infty}\psi_{k}^{(n)}\,\mathcal{L}^{-1}\\{e^{-b\sqrt{\frac{p}{A}}}p^{k+n/2-1}\\}\,.$
(B.73)

The inverse Laplace transform of an arbitrary term in the series,
$e^{-b\sqrt{\frac{p}{A}}}p^{k+n/2-1}$, can now be determined for all values of
$n$ and $k$ using Eq. C.109. Performing the said inverse we arrive at our
sought after asymptotic formula,

$C_{2}=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\psi_{k}^{(n)}\,g(-b/\sqrt{A},k+n/2-1,\tau)P_{n}(\cos(\phi))$
(B.74)

where

$\displaystyle g(x,m,\tau)=$
$\displaystyle\;\left\\{\begin{array}[]{ll}\tau^{-m}\left(\dfrac{\,{}_{1}F_{1}\left(m+1;\frac{1}{2};-\frac{x^{2}}{4\tau}\right)}{\tau\Gamma(-m)}-\dfrac{x\,\,_{1}F_{1}\left(m+\frac{3}{2};\frac{3}{2};-\frac{x^{2}}{4\tau}\right)}{\tau^{3/2}\Gamma\left(-m-\frac{1}{2}\right)}\right)&\;\textrm{if}\;\;m\leq-1,\\\
-\dfrac{x\,\tau^{-m-3/2}\,_{1}F_{1}\left(m+\frac{3}{2};\frac{3}{2};-\frac{x^{2}}{4\tau}\right)}{\Gamma\left(-m-\frac{1}{2}\right)}&\;\textrm{if}\;\;m=0,1,2,\ldots\,,\\\
\dfrac{\,\tau^{-m-1}\,_{1}F_{1}\left(m+1;\frac{1}{2};-\frac{x^{2}}{4\tau}\right)}{\Gamma(-m)}&\;\textrm{if}\;\;m=i+1/2,\;i=-1,0,1,2,\ldots\,,\end{array}\right.$
(B.78)

and ${}_{1}F_{1}$ is the confluent hypergeometric function of the first kind
[18].

The explicit forms of the $\psi_{k}^{(n)}$ are very large due to the
consecutive Cauchy products required for its computation. In the remainder of
this section we consider an approximation to Equation B.74, by considering
just the $k=0$ term of the inner summation. Hence, we get

$C_{2}\approx\sum_{n=0}^{\infty}\psi_{0}^{(n)}\,g(-b/\sqrt{A},n/2-1,\tau)P_{n}(\cos(\phi))\,.$
(B.79)

Using Eqs. B.52, B.55, B.56, B.61, B.65, B.71 & C.98, we obtain

$\psi_{0}^{(n)}=\frac{\sqrt{\textrm{p}}\,2^{-n}\,A^{1-\frac{n}{2}}\,X^{1+2\,n}\,\rho^{n}\,{\left(1+2\,n\right)}}{{\Gamma}\left(\frac{1}{2}+n\right)\,{\left(A\,X^{1+2\,n}+n\,{\left(A+A\,X^{1+2\,n}+X^{1+2\,n}-1\right)}\right)}}.$
(B.80)

In Section 3.2 we examine a better approximation of Eq. B.74 by considering
$\psi_{k}^{(n)}$ for $k=0,1,2$.

###  B.2 Large-time, asymptotic approximation for the concentration outside
the cell

We find the expansion for $\overline{C}_{1}$ as $p\to 0$ in a similar way to
the procedure followed in B.1. We write Eqs. A.42 & A.43 using new quantities,

$\overline{C}_{1}=\sum_{n=0}^{\infty}\overline{U}_{n}P_{n}(\cos(\phi)),\quad\overline{U}_{n}=\Big{(}\frac{(2n+1)}{p}e^{-b\sqrt{\frac{p}{A}}}\Big{)}\frac{\Omega_{3}}{\Omega_{2}}$
(B.81)

where

$\Omega_{3}=-(\sigma_{1}^{\prime}\omega_{3}-\sqrt{A}\sigma_{1}\omega_{3}^{\prime})\sigma_{5}\sigma_{4}+(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\sigma_{3}^{\prime})\sigma_{5}\omega_{4},\quad$
(B.82)

and where $\Omega_{2}$ is defined by Eq. B.46. To determine an expansion for
$\Omega_{3}$, we first substitute Eqs. B.48 & B.49 into the above expression
of $\Omega_{3}$ to get,

$\Omega_{3}=\frac{(-1)^{n+1}\textrm{p}}{2}((\sigma_{1}^{\prime}\beta_{3}-\sqrt{A}\sigma_{1}\beta_{3}^{\prime})\sigma_{4}-(\sigma_{1}^{\prime}\sigma_{3}-\sqrt{A}\sigma_{1}\,\sigma_{3})\beta_{4})\sigma_{5}\,.$
(B.83)

Then using the expansions defined by Eqs. B.51, B.53 & B.54, and evaluating
the product of these expansions using the quadruple Cauchy product (Eq.
C.105), we get

$\Omega_{3}=p^{n-1}\sum_{k=0}^{\infty}\alpha_{k}\,p^{k}$ (B.84)

where

$\displaystyle\alpha_{k}=\frac{(-1)^{n+1}\textrm{p}}{2}\times$
$\displaystyle\;\sum_{i=0}^{k}\sum_{j=0}^{i}\sum_{l=0}^{j}\sigma_{5,k-i}((\sigma_{1,i-j}^{\prime}\beta_{3,j-l}-\sqrt{A}\sigma_{1,i-j}\beta_{3,j-l}^{\prime})\sigma_{4,l}-(\sigma_{1,i-j}^{\prime}\sigma_{3,j-l}-\sqrt{A}\sigma_{1,i-j}\,\sigma_{3,j-l})\beta_{4,l}),$
(B.85)

where $\sigma_{l,k}$, $\sigma_{l,k}^{\prime}$, $\beta_{l,k}$ and
$\beta_{l,k}^{\prime}$ are defined by Eqs. B.52, B.55 & B.56. Substituting
Eqs. B.63 & B.84 into Eq. B.81, we have

$\overline{U}_{n}=(2n+1)p^{n/2-1}e^{-b\sqrt{\frac{p}{A}}}\Big{(}\sum_{k=0}^{\infty}\alpha_{k}\,p^{k}\Big{)}\Big{(}\sum_{k=0}^{\infty}\mu_{k}\,p^{k}\Big{)}$
(B.86)

where $\mu_{k}$, again appears, defined by Eq. B.63, the first four terms of
which are explicitly stated in Eq. B.65. Evaluating Eq. B.86 using the Cauchy
product,

$\overline{U}_{n}=p^{n/2-1}e^{-b\sqrt{\frac{p}{A}}}\sum_{k=0}^{\infty}\nu_{k}^{(n)}\,p_{k}$
(B.87)

where the coefficients

$\nu_{k}^{(n)}=(2n+1)\sum_{i=0}^{k}\alpha_{k-i}\,\mu_{i}\,,$ (B.88)

are independent of $p$. We are now in a position to determine the inverse
Laplace transform of $\overline{C}_{1}$ using Eq. B.87. From Eq. B.81, we have

$\mathcal{L}^{-1}\\{\overline{C}_{1}\\}=\sum_{n=0}^{\infty}\mathcal{L}^{-1}\\{\overline{U}_{n}\\}P_{n}(\cos(\phi)).$
(B.89)

where

$\mathcal{L}^{-1}\\{\overline{U}_{n}\\}=\sum_{k=0}^{\infty}\nu_{k}^{(n)}\,\mathcal{L}^{-1}\\{e^{-b\sqrt{\frac{p}{A}}}p^{k+n/2-1}\\}\,.$
(B.90)

The inverse Laplace transform of $e^{-b\sqrt{\frac{p}{A}}}p^{k+n/2-1}$ can be
determined for all values of $n$ and $k$ using Eq. C.109. Hence we have

$C_{1}=\sum_{n=0}^{\infty}\sum_{k=0}^{\infty}\nu_{k}^{(n)}\,g(-b/\sqrt{A},k+n/2-1,\tau)P_{n}(\cos(\phi))$
(B.91)

where $g(x,m,\tau)$ is again defined by Eq. B.78. The explicit forms of the
series coefficients, $\nu_{k}^{(n)}$, for arbitrary $k$ and $n$ are very large
due the consecutive Cauchy products required in the computation. To conclude
this section, however, we may consider a simple approximation to Equation B.91
by considering just the $k=0$ term for the inner summation. Hence, we get

$C_{1}\approx\sum_{n=0}^{\infty}\nu_{0}^{(n)}\,g(-b/\sqrt{A},n/2-1,\tau)P_{n}(\cos(\phi))\,.$
(B.92)

Using Eqs. B.52, B.55, B.56, B.61, B.65, B.85, B.88 & C.98, we find that

$\nu_{0}^{(n)}=\frac{\sqrt{\textrm{p}}\,2^{-1-n}\,X^{1+2\,n}\,{\left(n\,{\left(A+A\,\rho^{1+2\,n}+\rho^{1+2\,n}-1\right)}+A\,\rho^{1+2\,n}\right)}\,{\left(1+2\,n\right)}}{A^{n/2}\,\rho^{1+n}\,{\Gamma}\left(\frac{3}{2}+n\right)\,{\left(A\,X^{1+2\,n}+n\,{\left(A+A\,X^{1+2\,n}+X^{1+2\,n}-1\right)}\right)}}\,.$
(B.93)

In Section 3.2 we examine a better approximation of Eq. B.91 by considering
$\nu_{k}^{(n)}$ for $k=0,1,2$.

.

##  Appendix C Further Required Results

###  C.1 Modified spherical Bessel functions

The modified spherical Bessel functions of the first and second kind are
$i_{n}(z)$ ( $i_{-n}(z)$ is of negative order) and $k_{n}(z)$. We use the
definitions used by Olver et al. [23],

$i_{n}(z)=\sqrt{\frac{\textrm{p}}{2z}}I_{n+1/2}(z),\qquad
i_{-n}(z)=\sqrt{\frac{\textrm{p}}{2z}}I_{-n-1/2}(z),\qquad
k_{n}(z)=\sqrt{\frac{\textrm{p}}{2z}}K_{n+1/2}(z),$ (C.94)

where $I_{m}(z)$ and $K_{m}(z)$ are the modified Bessel functions of the first
and second kind of order $m$. We can relate $k_{n}(z)$ to $i_{n}(z)$ and
$i_{-n}(z)$ using the formula [23] (Eq. 10.47.11),

$k_{n}(z)=\frac{(-1)^{n+1}\textrm{p}}{2}(i_{n}(z)-i_{-n}(z)).$ (C.95)

The Wronskian of $i_{n}(z)$ and $k_{n}(z)$ is known to be [23] (Eq. 10.50.2),

$i_{n}^{\prime}(z)k_{n}(z)-i_{n}(z)k_{n}^{\prime}(z)=-\frac{\textrm{p}}{2z^{2}}\,.$
(C.96)

###  C.2 Series expansions

The series expansions of $i_{n}(z)$ and $i_{-n}(z)$ in [23] (Eq. 10.53.3 &
10.53.4) are here denoted

$i_{n}(z)=z^{n}\sum^{\infty}_{k=0}a(n,k)z^{2k},\qquad
i_{-n}(z)=z^{-n-1}\sum^{\infty}_{k=0}b(n,k)z^{2k}$ (C.97)

where

$a(n,k)=\frac{1}{2^{k}k!(2n+2k+1)!!},\quad
b(n,k)=\left\\{\begin{array}[]{ll}\dfrac{(-1)^{n+k}(2n-2k-1)!!}{2^{k}k!},&\quad\text{for}\quad
k\leq n,\\\ \dfrac{1}{2^{k}k!(2k-2n-1)!!},&\quad\text{for}\quad
k>n.\end{array}\right.$ (C.98)

Similarly, we have expressions of $i_{n}^{\prime}(z)$ and
$i_{-n}^{\prime}(z)$,

$\displaystyle i_{n}^{\prime}(z)=\frac{\partial(i_{n}(z))}{\partial z}$
$\displaystyle
i_{n}^{\prime}(z)=z^{n-1}\sum^{\infty}_{k=0}(2k+n)a(n,k)z^{2k},$ (C.99)
$\displaystyle i_{-n}^{\prime}(z)=\frac{\partial(i_{-n}(z))}{\partial z}$
$\displaystyle
i_{-n}^{\prime}(z)=z^{-n-2}\sum^{\infty}_{k=0}(2k-n-1)b(n,k)z^{2k},$ (C.100)

where the coefficients $a(n,k)$ and $b(n,k)$ are defined by Eq. C.98.

###  C.3 Multiple Cauchy products

We determine the product of three series. Let,

$P_{3}=\Big{(}\sum^{\infty}_{k=0}a_{k}\,z^{k}\Big{)}\Big{(}\sum^{\infty}_{k=0}b_{k}\,z^{k}\Big{)}\Big{(}\sum^{\infty}_{k=0}c_{k}\,z^{k}\Big{)}.$
(C.101)

Considering just the last two series,

$P_{3}=\Big{(}\sum^{\infty}_{k=0}a_{k}\,z^{k}\Big{)}\Big{(}\sum^{\infty}_{k=0}e_{k}\,z^{k}\Big{)},\quad
e_{k}=\sum^{k}_{j=0}b_{k-j}\,c_{j}\,,$ (C.102)

where $e_{k}$ was determined by the usual Cauchy product.Using this operation
again, we get

$P_{3}=\sum^{\infty}_{k=0}\Big{(}\sum^{k}_{i=0}a_{k-i}\,e_{i}\Big{)}\,z^{k}.$
(C.103)

Substituting Eq. C.102, we obtain the triple Cauchy product formula

$P_{3}=\sum^{\infty}_{k=0}\Big{(}\sum^{k}_{i=0}\sum^{i}_{j=0}a_{k-i}\,b_{i-j}\,c_{j}\Big{)}\,z^{k}\,.$
(C.104)

Similarly, for the case of a product of four series, we have the quadruple
Cauchy product formula,

$\Big{(}\sum^{\infty}_{k=0}a_{k}z^{k}\Big{)}\Big{(}\sum^{\infty}_{k=0}b_{k}z^{k}\Big{)}\Big{(}\sum^{\infty}_{k=0}c_{k}z^{k}\Big{)}\Big{(}\sum^{\infty}_{k=0}d_{k}z^{k}\Big{)}=\sum_{k=0}^{\infty}\Big{(}\sum_{i=0}^{k}\sum_{j=0}^{i}\sum_{l=0}^{j}a_{k-i}\,b_{i-j}\,c_{j-l}\,d_{l}\Big{)}\,z^{k}\,.$
(C.105)

###  C.4 Laplace inverse of expansions

Equation C.109 below is used to document the inverse Laplace transform of the
expansions developed in B. We have,

$\displaystyle\mathcal{L}^{-1}\\{e^{-x\sqrt{p}}p^{m}\\}=$
$\displaystyle\;\left\\{\begin{array}[]{ll}\tau^{-m}\left(\dfrac{\,{}_{1}F_{1}\left(m+1;\frac{1}{2};-\frac{x^{2}}{4\tau}\right)}{\tau\Gamma(-m)}-\dfrac{x\,_{1}F_{1}\left(m+\frac{3}{2};\frac{3}{2};-\frac{x^{2}}{4\tau}\right)}{\tau^{3/2}\Gamma\left(-m-\frac{1}{2}\right)}\right)&\;\textrm{if}\;\;m\leq-1,\\\
-\dfrac{x\,\tau^{-m-3/2}\,_{1}F_{1}\left(m+\frac{3}{2};\frac{3}{2};-\frac{x^{2}}{4\tau}\right)}{\Gamma\left(-m-\frac{1}{2}\right)}&\;\textrm{if}\;\;m=0,1,2,\ldots\,,\\\
\dfrac{\,\tau^{-m-1}\,_{1}F_{1}\left(m+1;\frac{1}{2};-\frac{x^{2}}{4\tau}\right)}{\Gamma(-m)}&\;\textrm{if}\;\;m=i+1/2,\;i=-1,0,1,2,\ldots\,,\end{array}\right.$
(C.109)

where

$_{1}F_{1}(a;b;z)=\sum_{k=0}^{\infty}\frac{\Gamma(a+k)\Gamma(b)}{\Gamma(b+k)\Gamma(a)\Gamma(k+1)}z^{k},$
(C.110)

$x>0$, $p$ is the Laplace parameter and ${}_{1}F_{1}$ is the confluent
hypergeometric function of the first kind [24]. Equation C.109 was determined
using the function `InverseLaplaceTransform` by Wolfram Mathematica [25]. The
expression condenses to a more comprehensible expressions for specific values
of $m$. For instance, $m=-1$

$\mathcal{L}^{-1}\\{e^{-x\sqrt{p}}p^{-1}\\}=\,_{1}F_{1}\left(0;\frac{1}{2};-\frac{x^{2}}{4\tau}\right)-\frac{x\,\,_{1}F_{1}\left(\frac{1}{2};\frac{3}{2};-\frac{x^{2}}{4\tau}\right)}{\sqrt{\tau\textrm{p}}}.$
(C.111)

Using formulas in [23] (Eq. 13.6.3 & 13.6.7) we observe

$\mathcal{L}^{-1}\\{e^{-x\sqrt{p}}p^{-1}\\}=\text{erfc}\left(\frac{x}{2\sqrt{\tau}}\right),$
(C.112)

which agrees with the expressions in [1, 21].