Hole-closing model reveals exponents for nonlinear degenerate diffusivity functions in cell biology

TL;DR

This paper investigates nonlinear degenerate diffusion models, specifically the Porous-Fisher equation, to better describe cell invasion and wound healing, demonstrating their ability to capture moving fronts and complex wound geometries.

Contribution

It introduces a hole-closing model using the Porous-Fisher equation, providing experimental and theoretical evidence supporting nonlinear degenerate diffusion in cell migration.

Findings

Nonlinear degenerate diffusion captures moving fronts in cell invasion.

Wound aspect ratio changes are consistent with degenerate diffusion models.

Linear diffusion fails to explain the elongation of wounds during closure.

Abstract

Continuum mathematical models for collective cell motion normally involve reaction-diffusion equations, such as the Fisher-KPP equation, with a linear diffusion term to describe cell motility and a logistic term to describe cell proliferation. While the Fisher-KPP equation and its generalisations are commonplace, a significant drawback for this family of models is that they are not able to capture the moving fronts that arise in cell invasion applications such as wound healing and tumour growth. An alternative, less common, approach is to include nonlinear degenerate diffusion in the models, such as in the Porous-Fisher equation, since solutions to the corresponding equations have compact support and therefore explicitly allow for moving fronts. We consider here a hole-closing problem for the Porous-Fisher equation whereby there is initially a simply connected region (the hole) with a nonzero population outside of the hole and a zero population inside. We outline how self-similar solutions (of the second kind) describe both circular and non-circular fronts in the hole-closing limit. Further, we present new experimental and theoretical evidence to support the use of nonlinear degenerate diffusion in models for collective cell motion. Our methodology involves setting up a 2D wound healing assay that has the geometry of a hole-closing problem, with cells initially seeded outside of a hole that closes as cells migrate and proliferate. For a particular class of fibroblast cells, the aspect ratio of an initially rectangular wound increases in time, so the wound becomes longer and thinner as it closes; our theoretical analysis shows that this behaviour is consistent with nonlinear degenerate diffusion but is not able to be captured with commonly used linear diffusion. This work is important because it provides a clear test for degenerate diffusion over linear diffusion in cell lines.

Figures

Equations

$$\frac{\partial u}{\partial t}=D\,\mbox{\boldmath$\nabla$}\cdot\left(\left(\frac{u}{K}\right)^{n}\mbox{\boldmath$\nabla$}u\right)+\lambda u\left(1-\frac{u}{K}\right).$$

$$\frac{\partial u}{\partial t}=D{\nabla }^{2}u+\lambda u\left(1-\frac{u}{K}\right),$$

$$u(\mathbf{x},0)=I(\mathbf{x}),\text{mbox}withI(\mathbf{x})=0\text{mbox}for\mathbf{x}\in \Omega (0),$$

$$\frac{\partial u}{\partial t}=\mbox{\boldmath$\nabla$}\cdot\left(u^{n}\mbox{\boldmath$\nabla$}u\right)+u\left(1-u\right)\quad\mbox{in}\quad{\bf x}\in\mathbb{R}^{2}\setminus\Omega(t),$$

$$u(r,0)=H(r-R(0)),$$

$$u=u^n\frac{\partial u}{\partial \nu }=0\text{mbox}on\partial \Omega (t),$$

$$\frac{\partial u}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}^n\frac{\partial u}{\partial r}\right)+u\left(1-u\right),r>R(t),$$

$$u=u^n\frac{\partial u}{\partial r}=0\text{mbox}onr=R(t),$$

$$u\sim (t_c-t{)}^{\alpha }U(\rho ),\mathrm{where}\rho =\frac{r}{(t_c-t{)}^{\beta }},$$

$$\frac{\partial u}{\partial t}=(t_c-t{)}^{\alpha -1}\left(-\alpha U+\beta \rho \frac{\mathrm{d}U}{\mathrm{d}\rho }\right),$$

$$\frac{\partial u}{\partial r}=(t_c-t{)}^{\alpha -\beta }\frac{\mathrm{d}U}{\mathrm{d}\rho },$$

$$\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}^n\frac{\partial u}{\partial r}\right)=(t_c-t{)}^{\alpha (n+1)-2\beta }\frac{1}{\rho }\frac{\mathrm{d}}{\mathrm{d}\rho }\left(\rho U^n\frac{\mathrm{d}U}{\mathrm{d}\rho }\right).$$

$$\frac{\partial u}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}^n\frac{\partial u}{\partial r}\right),r>R(t).$$

$$\frac{\partial u}{\partial t}\sim \frac{1}{r}\frac{\partial }{\partial r}\left(r{u}^n\frac{\partial u}{\partial r}\right)+u,r>R(t);$$

$$\tilde{u}={\mathrm{e}}^{(t_c-t)}u,\tilde{t}=\frac{1}{n}\left({\mathrm{e}}^{-n(t_c-t)}-1\right)$$

$$\frac{\partial \tilde{u}}{\partial \tilde{t}}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{\tilde{u}}^n\frac{\partial \tilde{u}}{\partial r}\right),r>R(\tilde{t}).$$

$$\alpha =\frac{2\beta -1}{n},$$

$$-\left(\frac{2\beta -1}{n}\right)U+\beta \rho \frac{\mathrm{d}U}{\mathrm{d}\rho }=\frac{1}{\rho }\frac{\mathrm{d}}{\mathrm{d}\rho }\left(\rho U^n\frac{\mathrm{d}U}{\mathrm{d}\rho }\right),\rho >\mu .$$

$$\frac{\partial u}{\partial t}=-\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}^n\frac{{\partial }^{3}u}{\partial r^3}\right),r>R(t),$$

$$-\left(\frac{2\beta -1}{n}\right)U+\beta \rho \frac{\mathrm{d}U}{\mathrm{d}\rho }\sim 0\text{mbox}as\rho \to \infty ,$$

$$U\sim a{\rho }^{(2\beta -1)/n\beta }\text{mbox}as\rho \to \infty .$$

$$U\to ϵU,\rho \to {ϵ}^{n/2}\rho ,$$

$$U\sim {\rho }^{(2\beta -1)/n\beta }\text{mbox}as\rho \to \infty .$$

$$U={\rho }_{\infty }^{(2\beta -1)/n\beta },\text{mbox}on\rho ={\rho }_{\infty },$$

$$\frac{\mathrm{d}U}{\mathrm{d}\rho }=\left(\frac{2\beta -1}{n\beta }\right){\rho }_{\infty }^{(2\beta -1-n\beta )/n\beta },\text{mbox}on\rho ={\rho }_{\infty }.$$

$$U=0\text{mbox}on\rho =\mu .$$

$$U^n\frac{\mathrm{d}U}{\mathrm{d}\rho }\to 0\text{mbox}as\rho \to {\mu }^{+}.$$

$$\beta \mu \frac{\mathrm{d}U}{\mathrm{d}\rho }\sim \frac{\mathrm{d}}{\mathrm{d}\rho }\left(U^n\frac{\mathrm{d}U}{\mathrm{d}\rho }\right)\text{mbox}as\rho \to {\mu }^{+}.$$

$$\beta \mu U\sim U^n\frac{\mathrm{d}U}{\mathrm{d}\rho }+\bar{C}\text{mbox}as\rho \to {\mu }^{+},$$

$$U\sim C(\rho -\mu {)}^{1/(n+1)}\text{mbox}as\rho \to {\mu }^{+}.$$

Full text

††journal: Physica D: Nonlinear Phenomena.

1 Introduction

The use of reaction-diffusion equations in mathematical biology is widespread [1, 2], especially for models of collective cell motion in applications of wound healing [3, 4, 5] and tumour invasion [6, 7, 8]. The most commonly used type of diffusion in these models is Fick’s first law, which gives rise to a linear diffusion term. However, a deficiency in this approach is that such governing equations do not allow for explicit descriptions of invading fronts. In an attempt to address this deficiency, some literature covering mathematical models of cell migration and proliferation has included studies of reaction-diffusion equations with nonlinear degenerate diffusion [9, 10, 11, 12, 13, 14, 15, 16]. These models allow for solutions with a well-defined moving boundary at the front of the invading cell population. An ongoing challenge in this area of research is to identify the most appropriate choice of nonlinear diffusion. We continue this work in the present paper by applying analytical and numerical techniques to study a hole-closing problem in the plane, focussing on the role of nonlinear degenerate diffusion and supporting our study with experimental results from an in vitro cell migration assay (that complements our recent study [17]).

In this paper we study the most simple nontrivial reaction diffusion equation with nonlinear degenerate diffusion that describes cell proliferation and migration, namely the Porous-Fisher equation

[TABLE]

This is a parabolic partial differential equation with nonlinear degenerate diffusion and a logistic growth term [18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. In cell biology, $u$ represents the density of a particular cell type, while ${\bf x}\in\mathbb{R}^{N}$, where $N=1$, 2 or 3. The diffusion term in (1) describes migration (or cell motility) of the cell population, with $n>0$ corresponding to a scenario in which cells are more likely to migrate if crowded, while the logistic growth component of (1) models cell proliferation with a carrying capacity $K$. The limiting case $n=0$ reduces (1) to the well-known the Fisher-Kolmogorov-Petrovski-Piskunov (Fisher-KPP) equation

[TABLE]

for which cell migration is due to random cell motility that is independent of cell density. As mentioned above, there are a plethora of studies in cell biology and ecology that use (2) or related models. In much of this vast literature, the Fisher-KPP equation has proved effective, and can successfully reproduce experimental behaviours; however, the choice of linear Fickian diffusion (or $n=0$ in (1)) is often also made because of its simplicity, not necessarily because of its biological relevance.

The main goal of the present study is to explore the role of nonlinear degenerate diffusion in a simple model of collective cell motion. As such, it is important to emphasise the key feature of the Porous-Fisher equation (1) with $n>0$ is that, unlike the Fisher-KPP equation (2) with linear diffusion, it allows solutions with compact support. This is a direct consequence of the term $(u/K)^{n}$, which vanishes in the limit $u\rightarrow 0^{+}$ (the diffusion is said to be degenerate in the sense that as $u\rightarrow 0^{+}$, the diffusion itself vanishes, and the equation changes from being of parabolic type to elliptic; to see this, write $v=u^{n}$ to give the eikonal equation \partial v/\partial t\sim|\mbox{\boldmath\nabla}v|^{2} in the limit $u\rightarrow 0^{+}$). Therefore, it is possible to use the Porous-Fisher equation (1) to model well-defined fronts of cell populations advancing on a region of zero population. On the other hand, even with initial conditions that have compact support, solutions to the Fisher-KPP equation (2) have $u>0$ for all ${\bf x}$ and $t>0$, meaning information is travelling infinitely fast, which strictly speaking is not biologically realistic.

The application we have in mind is a two-dimensional cell migration assay (${\bf x}\in\mathbb{R}^{2}$), with an invading population of cells moving over a substrate with a sharp front, ahead of which the cell population is essentially zero. In particular, we focus on the geometry in which there is initially a simply connected region devoid of cells, which we refer to as a hole or a wound. The resulting hole-closing problem is (1) subject to initial conditions

[TABLE]

where the simple closed curved $\partial\Omega(0)$ describes the initial shape of the hole (or wound). The challenge is to solve for $u$, but also to track the shape and speed of the boundary of the hole $\partial\Omega(t)$, especially in the hole-closing limit $t\rightarrow t_{c}^{-}$. We are motivated by our recent ([17]) and new experimental data from a wound healing (sticker) assay, as illustrated in Figure 1(a)-(b). In this particular example, the initial wound is circular in shape; however, a feature of our experimental design is that we are able to make the initial wounds any shape we choose. More generally, the study of hole-closing problems in cell biology has applications to a range of experimental wound healing scenarios, such as the circular wounds created by barrier assays [28, 29, 30], in a mammal’s ear [31] or cornea [32, 33] or a human skin equivalent construct [34], for example; similar mechanisms and geometry are involved in cell bridging experiments in tissue engineering [35]. Further motivation for our study arises from the lack of consensus when using (1) to model cell migration on the appropriate choice of the diffusion exponent $n$; we aim to shed light on this issue.

The outline of our paper is as follows. In the following section we set up our subsequent analysis by nondimensionaling the problem and presenting illustrative numerical solutions for the radially symmetric case of a circular hole (the numerical scheme is summarised in the Supplementary Material). In section 3, we argue that for times leading up to the hole-closing time $t_{c}$, the asymptotic behaviour of (1)-(3) in the neighbourhood of the point at which the hole closes ${\bf x}_{c}$ is equivalent to that for the Porous Medium equation (which is (1) with $\lambda=0$). For the radially symmetric case, we revisit the role of self-similarity of the second kind for solutions to the Porous Medium equation, reproducing some existing results by using an approach which is in some sense more transparent than that documented in the literature [36, 37, 38]. Numerical solutions are presented that illustrate how solutions to (1)-(3) with $k$-fold symmetry may or may not evolve to a circle in the hole-closing limit, depending on the diffusion exponent $n$, as per the analogous case of the Porous Medium equation [38, 39]. In section 4, we study solutions that are not self-similar in the hole-closing limit. In particular, we choose initial conditions for which the hole is initially rectangular and show how the hole becomes long and thin in the hole-closing limit, again following the behaviour of the Porous Medium equation [40]. By comparing with new experimental data from recently developed sticker assays (described in the Supplementary Material), we demonstrate how to estimate the exponent $n$ in (1). Finally, we close the paper in section 5 with a summary and a discussion about how our results provide support for the use of nonlinear degenerate diffusion in models for collective cell motion.

2 Governing equations and preliminary numerical results

2.1 Nondimensionalisation

As discussed in the Introduction, the model we focus on is the Porous-Fisher equation (1), which is the simplest model for collective cell motion in a two-dimensional assay that includes nonlinear degenerate diffusion. The only dependent variable is the cell density and, as a consequence, this model implicitly assumes that an excess of nutrients is available so that cell migration and cell proliferation are not affected by any lack of nutrient availability [41].

We nondimensionalise (1) by scaling cell density with respect to carrying capacity $K$ and choosing the representative length and time scales to be $\sqrt{D/\lambda}$ and $\lambda^{-1}$, respectively. As such, our dimensionless Porous-Fisher equation,

[TABLE]

does not involve any parameters apart from the diffusion exponent $n>0$. The dimensionless version of the initial condition (3) will, however, depend implicitly on the dimensional parameters. For example, as we discuss shortly, for a wound that is initially circular, the density $u=u(r,t)$ is a function of $r$ and $t$, and the relevant initial condition is

[TABLE]

where $H(r)$ is the Heaviside function and $R(0)$ is the initial radius of the wound. In this case, the dimensionless quantity $R(0)$ is the dimensional initial radius scaled by $\sqrt{D/\lambda}$.

We may formulate our hole-closing problem as a moving boundary problem by coupling (4) with the conditions

[TABLE]

where here $\partial u/\partial\nu$ is a normal derivative. The second condition in (6) enforces conservation of mass at the interface.

2.2 Numerical results for radially symmetric problem

We now provide numerical results for radially symmetric hole closing, which is governed by

[TABLE]

[TABLE]

subject to the initial condition (5). The numerical approach we use, based on applying a straight-forward finite-difference scheme on an uneven grid, is summarised in the Supplementary Material.

In Figure 1(c)-(f) we present density profiles for the four exponents $n=0.5$, $1$, $2$ and $3$, each with $R(0)=0.5$. In all four cases we see that the solutions have compact support so that $u>0$ for $r>R(t)$. Further, we note that the radius of the hole, $R(t)$, decreases in time until $R\rightarrow 0$ as $t\rightarrow t_{c}^{-}$, where $t_{c}$ is the hole-closing time. The dependence of $R(t)$ on $t$ is shown in Figure 1(g). We see here that the contact line $r=R(t)$ moves more slowly as the diffusion exponent $n$ increases, which is to be expected, since, for a fixed $u$, the nonlinear diffusion term $u^{n}$ decreases as $n$ increases.

One interesting observation is that the slope of the density profiles at the contact line $r=R(t)$ vanishes for $n<1$, is a nonzero constant for $n=1$, and is infinite for $n>1$ [24, 15]. Indeed, a simple leading order balance near the contact line for $R(t)>0$ suggests that $u\sim B(t)\,(r-R(t))^{1/n}$, which explains these qualitative behaviours.

Another observation about the propagation of the contact line is that $R(t)$ appears to follow a power-law $R(t)\sim\mu(t_{c}-t)^{\beta}$ in the hole-closing limit $t\rightarrow t_{c}^{-}$. This point is demonstrated in Figure 1(h), where the solid curves are numerical results corresponding to Figure 1(g). As these curves appear to approach a line on the log-log, a power-law behaviour is anticipated. The slope $\beta$ is dependent on $n$, as we discuss in detail in the following section.

3 Self-similar solutions

3.1 Similarity solutions for radially symmetric geometry

In the limit the hole closes, $t\rightarrow t_{c}^{-}$, we look for a similarity solution

[TABLE]

where at this stage the exponents $\alpha$ and $\beta$ are unknown. To proceed we require the partial derivatives

[TABLE]

[TABLE]

[TABLE]

We see that $\partial u/\partial t=\mathcal{O}((t_{c}-t)^{\alpha-1})$, $u=\mathcal{O}((t_{c}-t)^{\alpha})$ and $u^{2}=\mathcal{O}((t_{c}-t)^{2\alpha})$. Therefore, the source terms $u$ and $u^{2}$ in (7) do not contribute to leading order, and instead the self-similar behaviour is driven by nonlinear diffusion via the Porous Medium equation

[TABLE]

Furthermore, we even can keep the first source term $u$ so that

[TABLE]

in this case, the change of variables

[TABLE]

leads again to the Porous Medium equation

[TABLE]

We shall continue with (13) but simply note that our analysis of the Porous Medium equation holds exactly for (14).

By substituting (10) and (12) into (13), we find that

[TABLE]

[TABLE]

Here the similarity exponent $\beta$ cannot be determined by dimensional analysis of the governing equation (13) or by applying global conservation of mass; instead, it acts as an eigenvalue for a boundary-value problem associated with (18), which we discuss below in some detail. This self-similarity of the second kind makes the problem rather challenging.

3.2 Relationship to other studies

At this point it is worth making two comments on this formulation. First, the hole-closing problem for the Porous Medium equation (13) has been studied by a number of authors [36, 37, 38, 39] (including for the special case $n=3$ [42, 43], which corresponds to an inwardly-filling viscous gravity current). In all of these studies, the governing equation (13) is rewritten using the so-called pressure variable $v=u^{n}$ and the subsequent analysis involves similarity solutions of the form $v=(t_{c}-t)^{2\beta-1}V(\rho)$. Further changes of variables are required to complete the analysis, which to a certain extent has the effect of hiding the physical interpretation.

Second, the hole-closing problem for the fourth-order analogue of (13), the thin film equation

[TABLE]

has been studied recently by Zheng et al. [44, 45]. This analogue shares some features of the hole-closing problem for (13). In the spirit of [44], and in an attempt to provide a more transparent analysis than that which uses the pressure variable $v=u^{n}$, we shall concentrate on the dependent variable $u$ (and not $v$) and note some similarities and differences with [44] later.

3.3 Far-field conditions

To formulate the appropriate boundary-value problem associated with (18), we first derive the far-field condition as $\rho\rightarrow\infty$. As we do not want ${\partial u}/{\partial t}\rightarrow\infty$ as $t\rightarrow t_{c}^{-}$ for a fixed $r>0$ (in other words, our solution cannot blow up at any point other than $r=0$), we conclude from (10) and (17) that

[TABLE]

which implies

[TABLE]

This argument is used extensively in dealing with similarity solutions [46, 47]. The constant $a$ in (21) is arbitrary; however, the differential equation (18) is invariant under the transformation

[TABLE]

thus, without loss of generality, we can set $a=1$ and recover any solution we like by stretching $U$ and $\rho$ appropriately. Therefore we have as our far-field condition,

[TABLE]

Numerically, we will need to truncate $0<\rho<\infty$ to $0<\rho<\rho_{\infty}$, and then interpret (23) as two boundary conditions

[TABLE]

[TABLE]

Thus we can treat (18) with (24)-(25) as an initial-value problem starting at $\rho=\rho_{\infty}$.

3.4 Shooting method and near-field conditions

For a fixed $n$, we have a one-parameter family of initial-value problems, each for a different value of $\beta$. Our strategy is to interpret these via a shooting method, where we start at $\rho=\rho_{\infty}$ and shoot backwards until

[TABLE]

We wish to determine the appropriate value of $\beta$ for which (8) is satisfied, or alternatively,

[TABLE]

In addition to determining the appropriate value of $\beta$, the parameter $\mu$ must also be computed as part of the numerical solution. It will turn out that for each value of $n$, there are infinitely many pairs $(\mu,\beta)$ that satisfy (24)-(25) and (26), but only one pair $(\mu_{c},\beta_{c})$ which also satisfies (27). In order to proceed, we need to analyse possible behaviours of solutions to (18) as $U\rightarrow 0^{+}$.

Case 1, $\beta<\beta_{c}$. By considering (18) directly, we see that, provided $\mu\neq 0$, for (26) to hold we must have

[TABLE]

Integrating, we find

[TABLE]

where $\bar{C}$ is a constant. In the generic case where $\bar{C}\neq 0$, the two terms on the right-hand side balance, giving

[TABLE]

Following [44], we call this case generic touch-down.

Case 2, $\beta=\beta_{c}$. For the special case in (29) in which $\bar{C}=0$, the left-hand side must balance the first term on the right-hand side, giving

[TABLE]

Note in this case, (27) is satisfied, so this is the physically relevant solution we are after. We call this nongeneric touch-down [44].

Case 3, $\beta>\beta_{c}$. Here we now have $\mu=0$, which means that all the terms in (18) balance in the limit. By employing a power-law ansatz, we find that

[TABLE]

We call this touch-down at the origin [44].

In summary, the numerical task is to vary $\beta$ until the special borderline case 2 (nongeneric touch-down) is determined for $\beta=\beta_{c}$ and $\mu=\mu_{c}$. To ensure the borderline case is accurately identified, we utilise a high-order implicit finite-difference scheme to discretise (18), with local error tolerance set to near machine precision. The contrasting touch-down behaviours of cases 1 and 3 provide a convenient means of bracketing the critical value $\beta_{c}$, and hence convergence to this value is achieved iteratively through repeated bisection of the interval.

Some results of this task are presented in Figure 2. In Figure 2(a) the plots are for $n=0.5$. We can see clearly see the two representative profiles for $\beta<\beta_{c}$ (generic touch-down) have infinite slope at $\rho=\mu$; indeed, for this value of the diffusion exponent $n$, we have $U\sim C(\rho-\mu)^{2/3}$. On the other hand, the two representative profiles for $\beta>\beta_{c}$ (touch-down at the origin) are very flat in the limit (here $U\sim\rho^{4}/144$). The borderline case $\beta=\beta_{c}$ has $U$ scaling like $(\rho-\mu)^{2}$, which also has a zero slope, but does not approach the origin. The other three examples of $n$, shown in Figure 2(b)-(d), also show profiles for each of the three cases described above. We see that qualitatively the behaviour changes, depending on the value of the exponents in (30)-(32). In Figure 2(e)-(f) we show only the physically relevant similarity solutions with $\beta=\beta_{c}$. Again, the qualitative behaviour as these curves intersect the $\rho$-axis depends on the exponent in (31). In particular, we observe that the similarity solutions for $n=0.5$ and $0.75$ have zero slope at the contact line, while the solutions for $n=1.5$, $2$ and $3$ have infinite slope. The borderline case is $n=1$, where the slope is finite.

As a check on our similarity solutions, we present in Figure 3 profiles of $U$ versus $\rho$ (where we recall that $U=u(r,t)/(t_{c}-t)^{(2\beta-1)/n}$ and $\rho=r/(t_{c}-t)^{\beta}$) that are computed with our numerical scheme for various times. The figure includes plots for four values of $n$, namely $n=0.5$, $1$, $2$ and $3$. For each value of $n$, we see the initial condition is a scaled version of the Heaviside function (5), while as time increases, the numerical solutions of (7) (blue solid lines) approach the similarity solution (red dashed line) in the limit $t\rightarrow t_{c}^{-}$. This comparison between our numerical solution and our similarity solution provides confidence that our analysis in Section 3 is correct. Note we had to rescale our similarity solutions using (22) in order to match the point at which $U=0$ with the numerical solution.

We close this subsection by recalling the analogue problem for the thin film equation (19). As explained by Zheng et al. [44], the (circular) hole-closing problem for this fourth-order equation is similar to our problem (which is essentially for the Porous Medium equation (13)) in that there are self-similar solutions of the second kind. Further, our procedure for computing these similarity solutions is based on that presented in [44]. In particular, like in [44], we use a type of shooting method to explore numerical solutions for a range of parameter values and choose the relevant solution by ensuring the near-field limiting behaviour matches a physical constraint. The main difference between our study and that in [44] is that we have one free parameter $\beta$ and only three near-field options (which we call case 1, 2 and 3). On the other hand, the thin film equation is higher order and so Zheng et al.’s study is more complicated. They have a free parameter in the far-field condition as well as the similarity exponent $\beta$ and, consequently, they have more than three options for their near-field behaviour. Another further complication for the thin film model in [44] is that there needs to be a prewetting film which regularises a well-known singularity in stress at the moving contact line.

3.5 Similarity solutions with $k$-fold symmetry

For the Porous Medium equation, the stability of the radially symmetric similarity solutions (9) is known to depend on the diffusion exponent $n$. In particular, if we consider stability of the interface $R(t)=\mu(t_{c}-t)^{\beta}$ by adding a small perturbation $\gamma(t)\cos k\theta$, then for each $k\geq 3$, the perturbed solution is stable for $n>n_{k}$ and unstable for $n<n_{k}$, where $n_{k}$ is some borderline exponent that can be computed numerically [38, 39]. By stable, we mean that $\gamma/R\rightarrow 0$ as $t\rightarrow t_{c}^{-}$.

We postulate that the same type of stability holds for the full Porous Fisher equation (4), and now consider a relevant example for $k=4$. Suppose the initial condition is a square-shaped hole with a side of length unity. Here the solution has a 4-fold symmetry, so there are two options for the shape of the interface in the limit it closes. The first is for $n>n_{4}$, where we assume $n_{4}\approx 0.32$ [39]. For example, in Figure 4(a) we show the shape of the interface for a numerical solution with $n=0.5$. Since the radially symmetric similarity solution (9) is stable for $n=0.5$ to a 4-fold perturbation, we expect the interface for the full solution to approach a circle at extinction. While it is difficult to compute such two-dimensional solutions accurately near extinction, our numerical results in Figure 4(a) appears to show the hole becoming more circular as it closes. To support this idea, we have in Figure 4(b) plotted a kind of aspect ratio $\mathcal{A}_{\mathrm{d}}$, which is the ratio of the diagonal of the hole to the $x$-intercept. Again, it is difficult to tell given the scales involved, but it is not unreasonable to believe that $\mathcal{A}_{\mathrm{d}}$ is tending to unity at extinction, which should happen if it is becoming more circular.

Also shown in Figure 4 is an example for $n<n_{4}$, which demonstrates the second option for the shape of the hole as it closes. In Figure 4(c)-(d) the results are for $n=0.2$, for which we postulate the interface is unstable. Given the solution is rotationally symmetric (with 4-fold symmetry), this instability manifests itself by forcing the hole to approach a (noncircular) 4-fold symmetric shape at extinction, which is like a square with rounded corners. This type of evolution is suggested in Figure 4(c). Further, such a 4-fold symmetric shape will have a value of $\mathcal{A}_{\mathrm{d}}$ which is not unity. The time-dependence of $\mathcal{A}_{\mathrm{d}}$ in Figure 4(d) appears to support the idea that $\mathcal{A}_{\mathrm{d}}\not\to 1$ in the hole-closing limit.

We expect that qualitatively similar results could be presented for any $k$-fold symmetric initial condition, where $k\geq 3$, with the hole approaching a circle in shape for $n>n_{k}$ and a $k$-fold symmetric shape (which is like a regular $k$-sided polygon with rounded corners) for $n<n_{k}$. On the other hand, for $k=2$ we expect the similarity solution to unstable for all $n$, so that solutions with 2-fold symmetric initial conditions are no longer self-similar in the hole-closing limit [40, 38]; instead, the interface becomes oval in shape with an increasingly large aspect ratio which scales like $\mathcal{A}=\mathcal{O}((t_{c}-t)^{-1/2})$. We discuss this possibility further below in Section 4 where we focus on a rectangular-shaped initial hole.

All of these stability results are analogous to those for contracting bubbles in a Hele-Shaw cell for which the viscous fluid is of a power-law type [48, 49] or for which there is a competition between surface tension and a kinetic-type boundary condition [50, 51]. For these problems, there are also linear stability results that show a circular interface may be stable or unstable, depending on a parameter value, leading to the existence of noncircular self-similar solutions (that evolve to shapes which appear like $k$-sided polygons with rounded corners). Further, 2-fold symmetric perturbations are unstable for these Hele-Shaw problems, leading to interfaces that approach a slit in the hole-closing limit.

4 Rectangular shaped wounds

As just discussed, linear stability analysis of the radial similarity solutions shows that 2-fold symmetric perturbations grow in time [40], which suggests that an initially rectangular hole will become longer and thinner as it closes. We use this property to demonstrate how the diffusion exponent $n$ can be determined by fitting to aspect ratio data from a rectangular wound healing assay.

4.1 Sticker assays

We briefly summarise the sticker assays reported in [17]. These two-dimensional wound healing assays were performed using (NIH 3T3) fibroblast cells. A double-sided sticker was cut to a particular wound shape (circle, triangle or square) using a laser scribe, and then attached to a cell culture dish. The fibroblast cells were placed in the dish and incubated overnight. The sticker was removed to reveal the wound area of the required shape. As the cells migrated into the vacant space, images were taken at various discrete times and subsequently analysed using ImageJ [52] to determine the wound area at each time point (see Supplementary Material).

In [17], we simulated the sticker assays (with circular, triangular and square wounds), using a discrete random walk model on a hexagonal lattice incorporating crowding effects via an exclusion process [53, 54, 55] (whose continuum-limit description is the two-dimensional Fisher-KPP equation). We estimated the cell proliferation rate $\lambda=0.036$ /h and the carrying capacity $K=1.4\times 10^{-3}$ cells/$\mu$m2 by counting cells in sample regions of a corresponding proliferation assay and calibrating to the logistic growth model. The random walk model was then used to estimate the diffusion coefficient $D=1200\pm 260$ $\mu$m2/h which can be used for the Fisher-KPP equation.

4.2 New experimental results

We now report on new sticker assays performed for rectangular-shaped wounds using the same protocols as in [17]. For example, in Figure 5(a) we show experimental images for an initially rectangular wound whose aspect ratio is 2. These three representative images are taken for times $t=24$, 48 and 57 h. By approximating the wound boundary by a rectangle at each time step (see Supplementary Material), we are able to record the aspect ratio versus time for 3 replicates, and plot the result in Figure 5(d).

4.3 Fitting for the diffusion exponent $n$

To calibrate our experimental data with the Porous-Fisher equation (1), we first take the values $\lambda=0.036$ /h and $K=1.4\times 10^{-3}$ cells/$\mu$m2 from our previous study [17] described above. Then, for each fixed value of $n$ we choose to deal with, we fit for $D$ by comparing our numerical results of (7)-(8) with initial condition (5) with our previously obtained experimental data for the wound area with an initially circular wound [17]. We use a simple least-squares error to identify the most appropriate choice for $D$. Our results give estimates of $D$ and an interval of uncertainty in our estimates (Supplementary Material).

For each pair of $n$ and $D$ (together with $\lambda=0.036$ /h and $K=1.4\times 10^{-3}$ cells/$\mu$m2), we then solve the Porous-Fisher equation (1) numerically with a rectangular-shaped hole that represents the new experimental results shown in Figure 5(a). As time increases, the hole begins to close, as expected. Representative numerical results are shown in Figure 5(b)-(c) for $n=1$ and $2$ (note that the initial condition for these numerical solutions involves setting $u=1$ outside of the hole, truncated for numerical purposes at finite values of $x$ and $y$). Crucially, we plot the aspect ratio of the closing hole in Figure 5(e) for various values of $n$ and also include the experimental data in the same image. This strategy allows us to choose the value of $n$ which best matches the experimental data. A least-squares error between the experimental aspect ratio and the numerical results is shown in Figure 5(f). We see that, of the values of $n$ we have simulated, the choice $n=1$ appears to provide the closest match. Further details of the strategies employed in this subsection are provided in the Supplementary Material.

5 Discussion

In this paper we have studied various properties of solutions to the so-called hole-closing problem for the Porous-Fisher equation (1). We have summarised self-similar solutions of the second kind, which apply in the neighbourhood of the moving front in the limit the hole-closes. It turns out that the dynamics in this limit are governed by the Porous Medium equation and therefore many results carry over from previous studies of that equation [36, 37, 38, 39]. In contrast to those studies, we have formulated our analysis in terms of the original dependent variable and not the so-called pressure variable $v=u^{n}$; in this way our approach is more like the recent analogous study of the thin film equation [44]. Stability analysis [40] shows that $k$-fold (with $k\geq 3$) symmetric initial hole-shapes (like regular polygons) may become circular in the hole-closing limit, or may evolve to other non-circular shapes that keep the symmetry, all depending on the value of the exponent $n$. Further, 2-fold symmetric initial conditions (like rectangles) will become long and thin in the limit. Our numerical solutions confirm these predictions.

The present study is inspired by our recently published experimental data from a two-dimensional sticker assay with circular, square and triangular shaped wounds [17]. These experiments are modelled by hole-closing problems as there is an initially vacant wound area which is ultimately closed up as the cells migrate inwards and proliferate to occupy the initial wound space. In this paper we have presented new experimental results from sticker assays with rectangular-shaped wounds. We find the aspect ratio of the wounds increases in time in a way that agrees with our model using the Porous-Fisher equation (1). For various values of $n$, we are able to fit for the diffusivity $D$ by minimising the error between the numerically computed aspect ratio and the corresponding experimental data. Our results suggest that, for this cell line, a reasonable estimate for the exponent $n$ is $n=1$. This result is compatible with other studies of cell migration using the Porous-Fisher equation [10, 15]. Furthermore, this estimate for $n$ is consistent with our previous theoretical prediction [14], which suggests that for cells that themselves have an aspect ratio of $N$, the appropriate choice of $n$ is $n=N-1$. While NIH 3T3 fibroblast cells are not at all the same shape, many have an aspect ratio of roughly two.

To put these results in context, we recall how difficult it has been to identify an appropriate choice of the diffusion exponent $n$ in (1) when fitting with data from wound healing assays or experiments with traditional wound shapes. For example, for (approximately) one-dimensional fronts that arise from scratch assays with PC-3 prostate cancer cells, Jin et al. [10] calibrated the Porous-Fisher equation to the data for the examples $n=0.5$, $1$, $2$, $3$ and $4$. While it was concluded that the choice $n=1$ outperforms the others, the evidence was not straightforward as the corresponding estimates for the diffusion coefficient $D$ varied greatly over a range of initial conditions. In their study of circular wounds, Sherratt & Murray [12] briefly compared numerical simulations of the Porous-Fisher (1) with $n=4$ with experimental results from rabbit ears [31] and attempted to fit the data for $D$. While they observed the fit was not impressive, they did not attempt to vary $n$ but instead added other features of the model (to account for biochemical mediators). A further example is the study of Sengers et al. [13], who concluded that MG-63 bone cancer cells spread out radially in a way that was well represented by the Porous-Fisher with $n=1$, but again they did not attempt to fit the data with other values of the diffusion exponent. We conclude that our approach of using rectangular-shaped wounds in a sticker assay has the attractive feature of providing an additional means to fit for the diffusion exponent $n$ which has thus far been missing in the literature.

It is worth reflecting on how poorly the Fisher-KPP equation (2) performs at identifying certain properties of moving fronts at they evolve. As discussed above (and is well known), the Fisher-KPP equation does not allow solutions with compact support and so struggles to model scenarios with well-defined fronts. In particular, we are concerned with experiments where fronts of cell populations invade an empty space, such as in our two-dimensional wound healing assays. For this type of experiment, in order to apply the Fisher-KPP equation one must arbitrarily nominate a level set $u=\bar{u}$ as representing the moving front and then track properties of that level set as it evolves. For example, we have solved (2) numerically for a case with a $2\times 1$ rectangular shaped wound and plotted in Figure 6 the aspect ratio of a number of different level sets versus time. In all cases, the aspect ratios start at approximately $\mathcal{A}=2$ and then increase monotonically with time which is generally consistent with the experimental data (also included in this figure for reference). However, the actual values the aspect ratios take are very different for each level set and clearly there is no obvious choice as to which level set is most appropriate. As such, as a predictive tool, solutions to the Fisher-KPP equation (2) are not at all useful for describing this particular experimental property.

One possible feature of our experiments that is not included in the model (1) is chemotaxis, whereby cells produce a chemical signal, or chemoattractant, with concentration $g$, which can promote directional motion as the cells preferentially move up or down a gradient of $g$ [56]. The inclusion of chemotaxis in mathematical models for wound healing and tumour growth is commonplace [57, 4, 5, 58]. An extension of our model (1) which incorporates chemotaxis could be

[TABLE]

[TABLE]

where $\chi$ is the chemotactic sensitivity coefficient, $D_{g}$ is the diffusivity of the chemotactic chemical, $k_{1}$ is the rate at which cells produce the chemotactic chemical, and $k_{2}$ is the rate at which the chemotactic chemical undergoes natural decay [59, 60]. One significant challenge would be to obtain reasonable estimates for the four additional parameters $D_{g}$, $\chi$, $k_{1}$ and $k_{2}$ by calibrating the solution of the coupled system (33)-(34) with the experimental data. We have not pursued this approach for two main reasons. First, we have already been able to obtain a good match with the data using the simpler model (1) and we were able to use this calibrated model to draw conclusions about the use of degenerate diffusion for cell migration. Second, as is often the case with two-dimensional wound healing assays, we have not taken any measurements of concentrations of chemoattractants; without such measurements, there is obvious metholodogy for estimating the four additional parameters in the extended model (33)-(34). For this more complicated modelling to be useful, we suggest experimentalists make such measurements.

In summary, we have provided a range of evidence to support the use of nonlinear degenerate diffusion via the Porous-Fisher equation (1) for problems involving invading fronts. In particular, by comparing simulations with aspect ratio data taken from sticker assays with rectangular wounds, our methodology provides a clear test for degenerate diffusion over linear diffusion in two-dimensional cell migration experiments.

Acknowledgements

This work is supported by the Australian Research Council (DP140100249, DP170100474) and the Taiwan Ministry of Science and Technology (MOST 106-2313-B-002-031-MY3). WJ is supported by a QUT Vice Chancellor’s Research Fellowship. SWM acknowledges many useful discussions with David Frances and Sean McElwain. The authors are grateful for the computational resources and technical support provided by QUT’s High Performance Computing and Research Support group. Finally, the authors appreciate the supportive comments of the anonymous referees.