## Abstract

Reaction–diffusion waves in multiple spatial dimensions advance at a rate that strongly depends on the curvature of the wavefronts. These waves have important applications in many physical, ecological and biological systems. In this work, we analyse curvature dependences of travelling fronts in a single reaction–diffusion equation with general reaction term. We derive an exact, non-perturbative curvature dependence of the speed of travelling fronts that arises from transverse diffusion occurring parallel to the wavefront. Inward-propagating waves are characterized by three phases: an establishment phase dominated by initial and boundary conditions, a travelling-wave-like phase in which normal velocity matches standard results from singular perturbation theory and a dip-filling phase where the collision and interaction of fronts create additional curvature dependences to their progression rate. We analyse these behaviours and additional curvature dependences using a combination of asymptotic analyses and numerical simulations.

### 1. Introduction

Travelling fronts in reaction–diffusion models describe the progression in space of the transition between two states of an excitable system. These fronts represent the advance of a wave of excitation, where spatial locations in a rest state transition into an excited state [1,2]. Examples of such propagating waves abound in physical, biological and ecological systems. Particular examples of these waves are electrical activity along axonal membranes, chemical activity in chemical reactions, flame propagation in wildfires, biochemical activity in multicellular systems, biological tissue growth, biological species invasion, infection disease spread, social outbursts, as well as crystal growth and star formation in galaxies [1–10].

Wave propagation in these systems results from the combination of diffusion and autocatalytic production of the excitable element. In one spatial dimension, excitable media support the establishment of travelling waves with constant propagation speed $c$ in the long-time limit, under appropriate initial and boundary conditions [11]. In higher spatial dimensions, it is well known experimentally and theoretically that the propagation speed of travelling fronts is modified by the front’s curvature, which in general is a function of space and time [1,3,6,9,12]. Singular perturbation theories of excitable reaction–diffusion systems show that the normal velocity $v$ of travelling fronts is given by

In principle, equation (1.1) enables the determination of the spatial location of the travelling front as a function of time from an initial condition. It is referred to as the ‘eikonal equation for reaction–diffusion systems’ for this reason, by analogy to geometrical optics. Travelling front velocities are used in numerical marker particle methods as well as level-set methods applied to moving boundary problems [13–17]. Several results are known about the behaviour of fronts moving with curvature-dependent speeds. Curve-shortening flows and mean curvature flows for which normal velocity is proportional to curvature, as in equation (1.1), round off the initial front shape in such a way that the front becomes circular before shrinking into a point [13,18–21]. Such rounding of initial shape is commonly observed in experimental systems [3,7,9,22]. Many variants of mean curvature flows, such as length-preserving flows, area-preserving flows and Willmore flows are proposed as effective models of the evolution of fronts or interfaces [13,14,21,23–26]. There is recent interest in elucidating curvature-dependent velocities in moving fronts that represent the growth of biological tissues due to their application in tissue engineering [9,22,27–32]. In this context, the eikonal equation represents an understanding of a biological growth law based on migratory and proliferative cellular behaviours [33–37].

For the eikonal equation (1.1) to be a useful representation of the evolution of excitable systems, it is important that it represents the normal speed of travelling fronts accurately. Equation (1.1) is derived rigorously under mathematical assumptions on solution profiles. However, it is unknown how well these assumptions hold in practice. The plane-wave travelling speed $c$ is proportional to ${D}^{1/2}$ [1,11], so the curvature-dependent term in equation (1.1) is only relevant for curvatures of order $\mathrm{O}({D}^{-1/2})$. This corresponds to slowly moving highly curved fronts, such as where circular wavefronts collide and lead to an accelerating phase in species invasion fronts [1,6]. Such colliding circular fronts, however, can be expected to not conform to the assumptions used to derive equation (1.1) at some point. It is also unclear how boundary conditions may impact the eikonal equation for fronts located near boundaries of the domain.

In this work, we investigate the eikonal equation describing the normal speed of travelling fronts non-perturbatively in $D$. For simplicity, we consider a reaction–diffusion system consisting of a single species only, with linear diffusion, and an arbitrary reaction term. We perform numerical simulations of inward and outward wave propagation in domains of different shapes. This allows us to initiate fronts of the solution with various curvature, including shapes with positive and negative curvatures. We also investigate the influence of different boundary conditions imposed on the solution at the domain boundaries.

Our numerical simulations suggest that equation (1.1) is only valid in practice in a restricted time window, away from boundaries, and before circular fronts collide. We show numerically and analytically that there are several contributions to the curvature dependence of the normal speed of travelling fronts. The linear dependence upon curvature in equation (1.1) is an exact, non-perturbative contribution that originates from the transverse component of diffusion (diffusion parallel to the wavefront). Further dependences on curvature, including one of the same order in diffusivity, are exhibited explicitly using large-curvature asymptotic analyses and linear models. These further dependences arise from the normal component of diffusion and from autocatalytic production of the species in a spatially restricted environment.

### 2. Model description and theoretical developments

We consider the single-species reaction–diffusion model

A source term $F(u)$ of particular interest is the logistic growth term $F(u)=\lambda u(1-u)$ for a normalized density $u\in [0,1]$, in which case equation (2.1) is the two-dimensional analogue of the Fisher–Kolmogorov–Petrovsky–Piskunov (Fisher–KPP) model [11,39]. In this model, $u=0$ represents the unstable, excitable rest state, and $u=1$ represents the stable, excited state. In infinite space and under appropriate initial and boundary conditions, this model is well known to lead to travelling-wave solutions where rest states progressively transition into excited states [11,39,40]. However, we do not assume necessarily that the source term $F(u)$ leads to travelling-wave solutions, nor that it possesses rest states and excited states. We will state results valid for general source terms $F(u)$ and will illustrate properties of fronts of the solution to equation (2.1) with linear source terms, $F(u)=\lambda u$ and $F(u)=\lambda (1-u)$. These linear models do not support travelling-wave solutions. The time-dependent speed of their travelling fronts enables us to illustrate important properties of the eikonal equation. We also consider a bistable source term $F(u)=\lambda u(u-a)(1-u)$ with $0<a<1$, a common model for population growth subject to a strong Allee effect. The rest state $u=0$ and excited state $u=1$ in this model are both stable, yet this model still exhibits travelling-wave transitions from $u=0$ to $u=1$ or from $u=1$ to $u=0$ depending on $a$ and on the shape and size of initial conditions [4,41].

#### (a) Speed of travelling fronts

A travelling front of $u(\mathit{r},t)$ is defined to be the time-dependent set of points $\mathit{r}$ in the domain $\Omega $ such that $u(\mathit{r},t)={u}_{c}$ for a constant value ${u}_{c}$ (figure 1). A convenient way to estimate the velocity of travelling fronts is to calculate the velocity of contours of $u$ at each location of the domain. Each location $\mathit{r}$ of the domain has a density $u(\mathit{r},t)={u}_{c}$ for a certain allowable value ${u}_{c}$ of the solution. The velocity of a contour ${u}_{c}$ at $\mathit{r}$ can be determined as in level-set methods [13,14]. In short, given a parametrization $s\mapsto \mathit{r}(s,t)$ of a contour line at time $t$, the velocity of a point $\mathit{r}(s,t)$ of the contour at fixed $s$ is ${\mathit{r}}_{t}$. Differentiating $u(\mathit{r}(s,t),t))={u}_{c}$ with respect to $t$ and using the fact that the unit normal to level sets of $u$ in the decreasing direction of $u$ is $\mathit{n}=-\mathbf{\nabla}u/|\mathbf{\nabla}u|$, one obtains

#### (b) Exact curvature dependence due to transverse diffusion

Equation (2.2) enables us to extract an exact, non-perturbative contribution of the curvature of travelling fronts $u(\mathit{r},t)={u}_{c}$ to their progression rate $v$. In appendix A, we show that the Laplacian can be decomposed locally into a normal component and a transverse component with respect to the local geometry of the contour $u(\mathit{r},t)={u}_{c}$. The normal component of the Laplacian is ${u}_{nn}={\mathit{n}}^{\text{T}}H(u)\hspace{0.17em}\mathit{n}={\mathrm{\partial}}^{2}u/\mathrm{\partial}{n}^{2}$, where $H(u)=\mathbf{\nabla}\mathbf{\nabla}u$ is the Hessian matrix of $u$ and $\mathrm{\partial}/\mathrm{\partial}n=\mathit{n}\mathbf{\cdot}\mathbf{\nabla}$. The transverse component of the Laplacian, ${\mathrm{\nabla}}^{2}u-{u}_{nn}$, is always proportional to the mean curvature $\kappa =\mathbf{\nabla}\mathbf{\cdot}\mathit{n}$ of the travelling front (appendix A), i.e. ${\mathrm{\nabla}}^{2}u-{u}_{nn}=-\kappa |\mathbf{\nabla}u|$, so that

We note here that our sign conventions for the unit normal vector $\mathit{n}$ and curvature $\kappa =\mathbf{\nabla}\mathbf{\cdot}\mathit{n}$ are such that $\kappa <0$ at concave portions of the domain boundary and $\kappa >0$ at convex portions of the domain boundary. This sign convention and interpretation of convexity carries over to travelling fronts $u(\mathit{r},t)={u}_{c}$ inside the domain by interpreting regions where $u>{u}_{c}$ as occupied, and regions where $u<{u}_{c}$ as empty. With this convention, inward-moving circular fronts of the models $F(u)$ we consider have negative curvature, and outward-moving circular fronts have positive curvature.

The expression $v(\kappa )$ in equation (2.4) matches the asymptotic expression (1.1) from singular perturbation theory if $w$ is the speed $c$ of the travelling wave in the corresponding one-dimensional problem on an infinite domain, when such waves exist [1,4,6]. The relationship between $w$ and $c$ can be assessed by comparing equation (2.5) with the speed of travelling fronts in one dimension, given from equation (2.2) by

#### (c) Numerical simulations

Numerical simulations of equation (2.1) presented in this work are based on a simple explicit finite difference scheme using forward Euler time-stepping and centred differences in space (FTCS). Accuracy of numerical results is checked by refining computational grids until convergence is observed. Higher order time-stepping schemes are possible (such as total variation diminishing Runge-Kutta methods), but curvature effects are known to be more sensitive to spatial discretization than temporal discretization [14,36].

For general domain shapes $\Omega $, we discretize a square domain of size $L\times L$ large enough to contain $\Omega $ using a regular grid with constant space steps $\mathrm{\Delta}x=\mathrm{\Delta}y=h=L/(N-1)$. We define

##### (i) Circular domain

For circular domains $\Omega $ of radius *R* with radially symmetric solutions $u(r,t)$, where $r=|\mathit{r}|$, we implement a FTCS scheme by first expressing equation (2.1) in polar coordinates and assuming circular symmetry. The reaction–diffusion equation becomes

Dirichlet boundary conditions at the circular boundary $\mathrm{\partial}\Omega $ are implemented by setting ${u}_{N}^{k}={u}^{\ast}$. Neumann boundary conditions are implemented by setting ${u}_{N}^{k}={u}_{N-2}^{k}$ in equation (2.11).

##### (ii) Curvature and velocity

Numerical estimates of the mean curvature $\kappa $ are calculated from:

##### (iii) Domain boundaries

Figure 1 illustrates inward-moving fronts and outward-moving fronts on domains $\Omega $ specified by the following boundaries:

(i) | A circle of radius 1; | ||||

(ii) | A square of side length 2; | ||||

(iii) | A smooth cross shape parametrized in polar coordinates by the radial curve $$R(\theta )={\mathrm{sin}}^{4}(\theta )+{\mathrm{cos}}^{4}(\theta );$$ | ||||

(iv) | A petal with three branches (referred to below as three petal), parametrized in polar coordinates by the radial curve $$R(\theta )={R}_{0}+a\mathrm{cos}(n\theta ),\phantom{\rule{1em}{0ex}}\text{where}{R}_{0}=0.65,a=0.35,n=3.$$ |

### 3. Results and discussion

We start by examining numerical simulations of the two-dimensional analogue of the Fisher–KPP model with inward travelling fronts in the domains specified earlier. In this model, the reaction term represents logistic growth,

Figure 3 shows snapshots at three different times of numerical simulations of equation (2.1) with logistic source term (3.1) on the circular domain and Dirichlet boundary conditions with ${u}^{\ast}=1$. The solution exhibits travelling fronts progressing inward into the circular pore space (figure 3*a*). Contour lines $u={u}_{c}$ are shown every 0.1 increment of ${u}_{c}$. The contour line ${u}_{c}=0.5$, corresponding to locations with maximum growth rate, is shown in red.

Numerical estimates of front travelling speed and front curvature obtained at each (discretized) location $\mathit{r}$ of the domain are shown as data points $(v(\mathit{r},t),\kappa (\mathit{r},t))$ in figure 3*b*. These data points are coloured by the value ${u}_{c}$ of the corresponding contour. At all times, there is a well-defined linear correlation $v(\kappa )$ between front speed and curvature. At $t=3.0$, several data points with low density match the asymptotic relationship (1.1), shown in green, reasonably well. Numerical estimates of $\kappa $ and or $v$ at very low density at $t=3.0$ may be less accurate as they correspond to points near the centre of the domain, where spatial variations in $u$ are small. A value of $\kappa $ less than $-30$ corresponds to radii of curvature less than about three grid steps $h$. High-density data points are close to the boundary, and they deviate significantly from the relationship (1.1). At times $t=7.0$ and $t=15.0$, the relationship $v(\kappa )$ remains mostly linear, but the slope of the relationship is time dependent, in contrast to equation (1.1).

The results shown in figure 3 suggest that the asymptotic growth law (1.1), where the linear dependence of $v$ upon $\kappa $ has slope $-D$ and where $w$ in equation (2.4) is given by the one-dimensional travelling-wave speed $c$, equation (3.2), is well satisfied only during a limited time period, sufficiently far away from the domain boundary and before circular fronts collapse to a point. Not all points of the domain may exhibit the behaviour (1.1) at the same time, or at all. Equations (2.5)–(2.6) show that $w$ is well approximated by $c$ if the solution profile in the normal direction $\mathit{n}$ is close to the one-dimensional travelling-wave profile. To clarify which points of the domain satisfy this criterion in figure 3, we compare the solution profiles corresponding to figure 3 with the one-dimensional travelling-wave profile in figure 4*a*. This comparison confirms that the asymptotic growth law (1.1) in figure 3*b* is well satisfied when and where the solution profile in the normal direction matches the one-dimensional travelling-wave profile. Figure 4*b* also shows numerical results of inward progressing fronts of the two-dimensional Fisher–KPP model obtained with Neumann boundary conditions at $|\mathit{r}|=1$.

#### (a) Wave propagation phases

From these simulations, we can identify three distinct phases of wave propagation during inward motion:

(1) | An | ||||

(2) | A | ||||

(3) | A |

Numerical simulations of inward progressing fronts of the Fisher–KPP model in the square pore shape, cross pore shape and three-petal pore shape with Dirichlet boundary conditions are shown in figures 5, 6 and 7, respectively. In all pores, travelling fronts tend to smooth their initial shape and to become circular towards the centre of the domain. This is consistent with the evolution of interfaces governed by mean curvature flow, whereby normal velocity depends linearly on curvature [18,19,24]. Figures 5*b*, 6*b* and 7*b* show that there is a clear correlation between $v$ and $\kappa $ that becomes increasingly linear as time progresses. However, as mentioned earlier, this relationship is strongly time dependent. The asymptotic relationship (1.1) between $v$ and $\kappa $ is only well approximated during a limited period of time, away from the boundaries, and before fronts meet and interact in the centre. As for the circular pore shape, we can identify an establishment phase, a travelling-wave-like phase where solution profiles are similar to the one-dimensional travelling-wave profile (figure 8) and where $v(\kappa )$ matches the asymptotic relationship (1.1) well, and a dip-filling phase, characterized by a linear relationship $v(\kappa )$ but with a time-dependent slope.

The establishment phase results from complex interactions among boundary conditions, initial conditions, reaction and diffusion. In biological invasions, the establishment phase represents the survival and growth of a small initial population comprising few individuals, before the population expands out in the surrounding environment. This phase typically depends on many environmental and demographic conditions [42]. In our simulations, the initial population is distributed over the pore boundary. The establishment phase corresponds to the survival and the growth of the density $u$ at every location of the boundary.

When a well-defined travelling-wave phase exists, transitions between the different phases may be estimated qualitatively from the width of the one-dimensional travelling wavefront, i.e. the characteristic length over which the solution profile in the normal direction transitions from $u\approx 0$ to $u\approx 1$, and the one-dimensional travelling-wave speed $c$. The establishment phase transitions into the travelling-wave phase when the solution profile approximates the shape of the one-dimensional travelling wavefront over its characteristic width, and this solution profile detaches from the domain boundary. The travelling-wave phase transitions into the dip-filling phase when the width of the front begins to overlap with similar fronts travelling in opposite directions. In the Fisher–KPP model, the width of the travelling wave can be taken to be $4c$, the inverse of the steepness of the solution at $u={u}_{c}$ [40]. This width is directly related to the travelling wave speed $c$.

By decreasing diffusivity, or equivalently, by increasing domain size [9], the travelling-wave-like phase lasts longer, and the match between $v(\kappa )$ and the asymptotic relationship (1.1) improves (figure 9). By contrast, if diffusivity is too large or domain size too small, the solution may never exhibit a well-defined travelling-wave-like phase. In figures 3 and 9, the one-dimensional travelling wavefront has the same width and speed. However, establishing a profile similar to that of the one-dimensional travelling wave takes longer in the larger circular pore and results in a longer establishment phase (figure 9*b*).

The normal velocity $v$ is not always a sole function of curvature. The spread of points $(\kappa ,v)$ in figures 6*b* and 7*b* indicates that other details of the solution may influence the normal velocity. The curves formed by the spread of points in these figures are due to only calculating $v$ and $\kappa $ at discretized positions in the domain. In continuous space, the spread of points would fill a whole region, similar to what is seen at low curvature (more points of the computational grid have small curvature than high curvature).

In figure 7, the front $u={u}_{c}$ splits into multiple disconnected fronts ($t=2.8$). Dip-filling occurs independently in the middle of the petals and in the centre of the domain. These separate dip-filling locations correspond to local maxima of the closest distance to the domain boundary [13].

Outward-moving fronts clearly do not exhibit a dip-filling phase (figure 1). The long-time solution $u(\mathit{r},t)$ develops a profile in the normal direction that gradually converges to that of the one-dimensional infinite space solution [12,42]. The asymptotic relationship (1.1) is well satisfied sufficiently far from the boundaries, but the correction term due to curvature is small and tends to zero as $t\to \mathrm{\infty}$. Indeed, as the solution expands out away from the origin, the curvature of the front decreases. In the following, we focus on analysing inward-moving fronts and the long-term behaviour of the solution in the dip-filling phase, since in this situation, the curvature of inward-moving fronts increases with time and induces further dependences to the front speed, as suggested by our numerical simulations.

#### (b) Normal diffusion in the dip-filling phase

Figures 3–7 suggest that the dip-filling phase is similar across many domain shapes and leads in the long term to a well-defined linear dependence $v(\kappa )$ of front speed upon curvature, although with a different slope compared with the asymptotic perturbation result (1.1). This similarity across domain shapes is due to the fact that the travelling fronts always evolve into shrinking circles in these examples, an evolution likely initially due to the explicit curvature dependence of normal velocity in equation (2.4) [18,19]. As time progresses, further curvature dependences of $v(\kappa )$ arise that reinforce this tendency.

To analyse these additional contributions to $v(\kappa )$ in the dip-filling phase, we now assume that the solution is approximately radially symmetric near the centre of the domain $r\to 0$. Since $\mathit{n}$ points in the opposite direction to the radial axis, $\mathrm{\partial}/\mathrm{\partial}n=\mathit{n}\mathbf{\cdot}\mathbf{\nabla}=-\mathrm{\partial}/\mathrm{\partial}r$. The contribution to $v(\kappa )$ due to diffusion in the normal direction captured by the term $D{u}_{nn}/|\mathbf{\nabla}u|$ in equation (2.5) becomes, using equations (2.9) and Bernoulli–L’Hôpital’s rule (2.10):

Importantly, the contribution to $v(\kappa )$ due to the normal component of diffusion is of the same order as the contribution due to the transverse component of diffusion obtained usually by singular perturbation analyses. These singular perturbation analyses thus do not describe the velocity of inward-moving circular fronts close to where they collide [1]; these developments are only valid away from boundaries and away from colliding fronts.

#### (c) Reaction-rate dependence in the dip-filling phase

The contribution due to the reaction term in equation (3.3) also possesses a dominant linear dependence upon curvature. Indeed, proceeding similarly to §3b and using ${u}_{r}/r\sim {u}_{rr}$ as $r\to 0$ and ${u}_{n}=-{u}_{r}$, we have

#### (d) Zero diffusion limit

Equations (3.3) and (3.5) show that the dip-filling phase can be highly influenced by reaction terms since fronts of the solution may progress in space due to $u$ increasing locally (figure 2). However, the reaction rate-dependent contribution to $v(\kappa )$ on the right-hand side of equation (3.5) is still coupled to diffusive effects. So, here, we investigate two remaining important points: (i) the influence of diffusion in the reaction rate-dependent contribution to $v(\kappa )$; and (ii) how far from the origin $r\to 0$ we can expect the asymptotic behaviour of $v(\kappa )$ in equation (3.5) to hold in practice.

To obtain insights into the velocity of fronts due to reaction terms only, we consider the zero diffusion limit $D\to 0$ in the dip-filling phase. We show that in this regime, the eikonal equation $v(\kappa )$ becomes time independent for any choice of $F(u)$. In other words, diffusive effects are entirely responsible for the time dependence of the slope of $v(\kappa )$ observed at long times in figures 3 and 5–7, and captured in equation (3.5). First, let us show that the velocity field $v(\mathit{r},t)$ of all moving fronts is independent of time if $D=0$. From equation (2.2) with $D=0$, we have

Figure 11 shows numerical simulations in which the two-dimensional Fisher–KPP model is evolved in a circular pore until time $t=9.0$, at which point diffusion is set to zero. From this point onwards, the solution continues to evolve due to logistic growth only (figure 11*a*). We observe that the relationship $v(\kappa )$ in the dip-filling phase is linear with a slope that becomes steeper with time until $t=9.0$. After this point, the relationship $v(\kappa )$ no longer evolves in time and is approximately given by $v(\kappa )=-\gamma \kappa $, where $\gamma >0$ is a constant (figure 11*b*).

The proportionality constant $\gamma $ is directly related to the profile of the solution at $t=9.0$ and its subsequent evolution at times $t>9.0$, as follows. Let ${t}_{0}$ be the time at which $D$ is set to zero. When $D=0$ and with radial symmetry (dip-filling phase), the solution profile at time $t>{t}_{0}$ is given by a function $U(r,t)$, where $r\ge 0$ is the radial coordinate. From equation (2.2) and assuming $v(\kappa )=-\gamma \kappa =\gamma /r$, we have:

The value of $\gamma $ that best fits the numerical solution profile at $t=9.0$ around $r=0,$ as shown in figure 11*a*, is $\gamma \approx 0.047$. This value of $\gamma $ found from the solution profile predicts the slope of the growth law $v(\kappa )$ very well, see figure 11*b*. The solution in equation (3.9) differs from the numerical solution at locations where fronts have small curvature ($|x|\gtrsim 0.2$, corresponding to $|\kappa |\lesssim 5$, see figure 11*a*). At these curvatures, numerical estimates of $v(\kappa )$ in figure 11*b* deviate from the linear relationship $v(\kappa )=-\gamma \kappa $. Figure 11*a* also shows that once diffusion is turned off, the time evolution of the solution is very well represented by equation (3.9).

The correspondence shown in equation (3.7) between the eikonal equation $v(\kappa )$ observed at a fixed time and solution profiles in the radial direction around the origin may be generalized to time-dependent and nonlinear relationships $v(\kappa )$ by integrating, from equation (2.4),

In summary, the dip-filling phase is characterized by travelling front velocities $v(\kappa )$ that are affected by curvature due to several contributions: (i) the transverse component of diffusion, always exactly contributing a term $-D\kappa $; (ii) the normal component of diffusion, contributing an extra term $-D\kappa $ at large curvature $\kappa \to -\mathrm{\infty}$; and (iii) reaction terms, which provide additional contributions to the normal speed that also becomes linear in $\kappa $ at large curvature $\kappa \to -\mathrm{\infty}$, with a slope whose time dependence is entirely due to diffusive effects.

#### (e) Linearized models

To understand in more detail how solution profiles converge in the dip-filling phase to profiles that precisely give rise to linear contributions in the growth law $v(\kappa )$, and to exhibit explicit time dependences of the slope of these linear contributions in the long-time limit, we now consider two relevant linear models: Skellam’s model, where $F(u)=\lambda u$, and a saturated growth model, where $F(u)=\lambda (1-u)$. Skellam’s model corresponds to linearizing the logistic reaction rate about the unstable rest state $u=0$ and is thus expected to describe the two-dimensional analogue of the Fisher–KPP model at early times. The saturated growth model corresponds to linearizing the logistic reaction rate about the stable excited state $u=1$ reached by the solution of the Fisher–KPP model in the long-time limit.

Since our focus is the dip-filling phase, where solutions develop circular fronts in the long term, we assume radially symmetric solutions and solve equation (2.1) for a radially symmetric solution $u(r,t)$ with

To solve the equivalent problem with the saturated growth reaction term

We can further elucidate the contribution to $v(\kappa )$ due to the reaction term, $F(u)/{u}_{r}$, in the long-time limit based on equations (3.11)–(3.14). In Skellam’s model with Dirichlet boundary conditions, we obtain

In the saturated growth model, by retaining the next-order term in the long-time limit, we obtain:

We see that in Skellam’s model, the slope of the linear relationship $v(\kappa )$ continues to increase with time (in absolute value). By contrast, in the saturated growth model, the slope of the linear relationship $v(\kappa )$ decreases with time (in absolute value) and converges to

#### (f) Strong Allee effect

Finally, we consider an example of bistable reaction term where both the rest state $u=0$ and excited state $u=1$ are stable steady states. A classic example is provided by population growth with strong Allee effect, modelled by $F(u)=\lambda u(u-a)(1-u)$ with $0<a<1$. This model supports travelling waves with speed

Our numerical observation that the linear relationship $v(\kappa )$ of the two-dimensional analogue of the Fisher–KPP model and of the strong Allee effect converge to that of the saturated growth model in the long time is not trivial, even if the reaction terms in these models becomes similar in this limit. From equation (3.5), the slope of the relationship $v(\kappa )$ depends on ${u}_{rr}(0,t)$. This quantity could in principle depend on how the solution profile approaches the steady state from the initial condition, where the reaction term has different contributions in the three models.

### 4. Conclusion and future work

Reaction–diffusion waves propagating in excitable systems in multiple spatial dimensions are strongly influenced by travelling front curvature. The eikonal equation for these parabolic systems of differential equations is an expression that determines the local normal velocity of travelling fronts. In principle, a complete description of the eikonal equation at each location of the domain and each time enables one to fully evolve the solution in time, much like the method of characteristics for hyperbolic conservation laws [44]. The eikonal equation thus provides valuable physical insight into the behaviour of reaction–diffusion systems, as well as mathematical insights from known results from geometric flows that evolve moving interfaces based on velocity laws [19,26]. In a tissue engineering context, the eikonal equation provides a biological growth law that describes how bulk tissue production and diffusive redistribution of tissue material lead to differential progression rates of tissue interfaces during tissue growth or invasion [9,32,33].

In the present work, we analysed the eikonal equation for a single reaction–diffusion equation with arbitrary reaction term. We showed that the contribution $-D\kappa $ to the eikonal equation derived from singular perturbation theory in the low-diffusion limit [1,4,6] is in fact an exact, non-perturbative contribution originating from diffusion occurring transversally, i.e. along the wavefront. For fronts moving inward into an empty domain, we identify three phases of the solution: (i) an establishment phase, strongly influenced by conditions on the domain boundary; (ii) a travelling-wave-like phase, well described by the eikonal equation $v=c-D\kappa $ predicted by singular perturbation theory in the low-diffusion limit; and (iii) a dip-filling phase, in which the eikonal equation possesses further curvature dependences arising from the collision and interaction of inward-moving fronts. At large curvature, the normal velocity of colliding circular fronts is proportional to curvature, with an additional contribution $-D\kappa $ originating from diffusion occurring normally to wavefronts, and a time-dependent contribution originating from non-conservative changes to the solution described by the reaction term. The latter contribution becomes time independent if $D=0$ for any choice of the reaction term $F(u)$. These behaviours are confirmed by the explicit long-time solutions of Skellam’s model ($F(u)=\lambda u$), the saturated growth model ($F(u)=\lambda (1-u)$), as well as the strong Allee source term ($F(u)=\lambda u(u-a)(1-u)$), which is a common model for bistable populations [4,41,45].

Our numerical simulations suggest that for a broad range of domains and reaction terms $F(u)$, inward-moving fronts that meet and interact evolve into colliding circular fronts. This behaviour is expected to be due initially to the exact linear dependence of $v$ upon $\kappa $ that arises from transverse diffusion, and the fact that mean curvature flows evolve interfaces into shrinking circles [18–21]. However, this behaviour likely depends on the strength of this curvature-dependent term relative to other contributions to the normal velocity. There are situations in reaction–diffusion systems with degenerate nonlinear diffusion where inward-moving fronts do not round off as circles before disappearing [30], in which case other expressions for the eikonal equation may hold. Excitable systems comprising several reaction–diffusion equations can present more complex wave patterns, such as spiral waves and wave trains [11]. Calcium waves and other reaction–diffusion waves may also evolve on the surface of curved domains [46,47]. In these cases, the curvature of the surface also influences travelling-wave speeds [2,48]. It is unclear how our results apply to these situations. Generalizing the methods presented in this article to more general diffusion operators including the Laplace–Beltrami operator could be the subject of future investigations. Another situation of interest is the investigation of curvature-dependent front speeds in metastable systems with wave-pinning, such as the Allen–Cahn equation with $F(u)=u(1-{u}^{2})$ [49], where the dip-filling phase is absent. Finally, generalizations of our analyses to heterogeneous excitable media [50] and to Fisher–Stefan models representing reaction–diffusion systems coupled with moving boundary conditions [51,52] are also of interest, in both situations supporting travelling-wave solutions as well as dip-filling situations.

### Data accessibility

Key algorithms used to generate results are freely available on Github at https://github.com/prbuen/Buenzli2022_CurvatureReacDiff.

### Authors' contributions

P.B.: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, software, validation, visualization, writing—original draft and writing—review and editing; M.J.S.: conceptualization, investigation, methodology and writing—review and editing.

Both authors gave final approval for publication and agreed to be held accountable for the work performed therein.

### Conflict of interest declaration

We declare we have no competing interests.

### Funding

This research was supported by the Australian Research Council (DP180101797 and DP190102545).

## Appendix A. Normal and transverse components of the Laplacian

In this appendix, we decompose the Laplacian into a component normal to level sets of $u$ and transverse components tangent to level sets of $u$ and show that the transverse components are proportional to the mean curvature $\kappa =\mathbf{\nabla}\mathbf{\cdot}\mathit{n}$. The mean curvature here is defined as the sum of principal curvatures, rather than their average [13,24], and $\mathit{n}=-\mathbf{\nabla}u/|\mathbf{\nabla}u|$ is the unit normal vector to level sets of $u$ pointing towards decreasing values of $u$. We start by expanding the divergence of $\mathit{n}$ using standard rules of differentiation:

### References

- 1.
Tyson JJ, Keener JP . 1988 Singular perturbation theory of traveling waves in excitable media (a review).**Physica D**, 327-361. (doi:10.1016/0167-2789(88)90062-0) Crossref, Web of Science, Google Scholar**32** - 2.
Grindrod P, Lewis MA, Murray JD . 1991 A geometrical approach to wave-type solutions of excitable reaction–diffusion systems.**Proc. R. Soc. Lond. A**, 151-164. (doi:10.1098/rspa.1991.0040) Link, Web of Science, Google Scholar**433** - 3.
Foerster P, Müller SC, Hess B . 1988 Curvature and propagation velocity of chemical waves.**Science**, 685-687. (doi:10.1126/science.241.4866.685) Crossref, PubMed, Web of Science, Google Scholar**241** - 4.
Lewis MA, Kareiva P . 1993 Allee dynamics and the spread of invading organisms.**Theor. Pop. Biol.**, 141-158. (doi:10.1006/tpbi.1993.1007) Crossref, Web of Science, Google Scholar**43** - 5.
Echekki T, Chen JH . 1996 Unsteady strain rate and curvature effects in turbulent premixed methane-air flames.**Combust. Flame**, 184-202. (doi:10.1016/0010-2180(96)00011-9) Crossref, Web of Science, Google Scholar**106** - 6.
Volpert V, Petrovskii S . 2009 Reaction–diffusion waves in biology.**Phys. Life Rev.**, 267-310. (doi:10.1016/j.plrev.2009.10.002) Crossref, PubMed, Web of Science, Google Scholar**6** - 7.
Hilton JE, Miller C, Sharples JJ, Sullivan AL . 2016 Curvature effects in the dynamic propagation of wildfires.**Int. J. Wildland Fire**, 1238-1251. (doi:10.1071/WF16070) Crossref, Web of Science, Google Scholar**25** - 8.
Dave HL, Chadhury S . 2020 Evolution of local flame displacement speedsin turbulence.**J. Fluid Mech.**, A46. (doi:10.1017/jfm.2019.896) Crossref, Web of Science, Google Scholar**884** - 9.
Buenzli PR, Lanaro M, Wong CS, McLaughlin MP, Allenby MC, Woodruff MA, Simpson MJ . 2020 Cell proliferation and migration explain pore bridging dynamics in 3D printed scaffolds of different pore size.**Acta Biomater.**, 285-295. (doi:10.1016/j.actbio.2020.07.010) Crossref, PubMed, Web of Science, Google Scholar**114** - 10.
Bakhshi M, Ghazaryan A, Manukian V, Rodgriguez N . 2021 Traveling wave solutions in a model for social outbursts in a tension-inhibitive regime.**Stud. Appl. Math.**, 650-674. (doi:10.1111/sapm.12394) Crossref, Web of Science, Google Scholar**147** - 11.
Murray JD . 2004**Mathematical biology: II Spatial models and biomedical applications**. New York, NY: Springer. Google Scholar - 12.
Witelski TP, Ono K, Kaper TJ . 2000 On axisymmetric traveling waves and radial solutions of semi-linear elliptic equations.**Natural Resource Modeling**, 339-388. (doi:10.1111/j.1939-7445.2000.tb00039.x) Crossref, Google Scholar**13** - 13.
Sethian J . 1999**Level set methods and fast marching methods**. Cambridge, UK: Cambridge University Press. Google Scholar - 14.
Osher S, Fedkiw R . 2003**Level set methods and dynamic implicit surfaces**. New York, NY: Springer. Crossref, Google Scholar - 15.
García-Reimbert C, Minzoni AA, Ondarza R, Herrera JJE . 2002 Simple diffraction processes for Nagumo- and Fisher-type diffusion fronts.**Stud. Appl. Math.**, 367-389. (doi:10.1111/1467-9590.1074191) Crossref, Web of Science, Google Scholar**107** - 16.
Leung S, Zhao H . 2009 A grid based particle method for moving interface problems.**J. Comp. Phys.**, 2993-3024. (doi:10.1016/j.jcp.2009.01.005) Crossref, Web of Science, Google Scholar**228** - 17.
Hon SY, Leung S, Zhao H . 2014 A cell based particle method for modeling dynamic interfaces.**J. Comp. Phys.**, 279-306. (doi:10.1016/j.jcp.2014.04.032) Crossref, Web of Science, Google Scholar**272** - 18.
Brakke KA . 1978**The motion of a surface by its mean curvature**. Princeton, NJ: Princeton University Press. Google Scholar - 19.
Grayson MA . 1987 The heat equation shrinks embedded plane curves to round points.**J. Differ. Geom.**, 285-314. (doi:10.4310/jdg/1214441371) Crossref, Web of Science, Google Scholar**26** - 20.
Chou K-S, Zhu X-P . 1999 A convexity theorem for a class of anisotropic flows of plane curves.**Indiana University Math. J.**, 139-154. (doi:10.1512/iumj.1999.48.1273) Crossref, Web of Science, Google Scholar**48** - 21.
Liu Z . 2013 Mean curvature flow with a constant forcing.**Nonlinear Differ. Equ. Appl.**, 621-649. (doi:10.1007/s00030-012-0171-4) Crossref, Web of Science, Google Scholar**20** - 22.
Rumpler M, Woesz A, Dunlop JWC, van Dongen JT, Fratzl P . 2008 The effect of geometry on three-dimensional tissue growth.**J. R. Soc. Interface**, 1173-1180. (doi:10.1098/rsif.2008.0064) Link, Web of Science, Google Scholar**5** - 23.
Chen Y-G, Giga Y, Goto S . 1991 Uniqueness and existence of viscosity solutions of generalized mean curvature flow equations.**J. Differ. Geom.**, 749-786. (doi:10.4310/jdg/1214446564) Crossref, Web of Science, Google Scholar**33** - 24.
- 25.
Dallaston MC, McCue SW . 2016 A curve shortening flow rule for closed embedded plane curves with a prescribed rate of change in enclosed area.**Proc. R. Soc. A**, 20150629. (doi:10.1098/rspa.2015.0629) Link, Google Scholar**472** - 26.
Andrews B, Chow B, Guenther C, Langford M . 2020**Extrinsic geometric flows**. Providence, RI: American Mathematical Society. Crossref, Google Scholar - 27.
Bidan CM, Kommareddy KP, Rumpler M, Kollmannsberger P, Fratzl P, Dunlop JWC . 2012 Geometry as a factor for tissue growth: toward shape optimization of tissue engineering scaffolds.**Adv. Healthcare Mater.**, 186-194. (doi:10.1002/adhm.201200159) Crossref, PubMed, Web of Science, Google Scholar**2** - 28.
Guyot Y, Luyten FP, Schrooten J, Papantoniou I, Geris L . 2015 A three-dimensional computational fluid dynamics model of shear stress distribution during neotissue growth in a perfusion bioreactor.**Biotech. Bioeng.**, 2591-2600. (doi:10.1002/bit.25672) Crossref, PubMed, Web of Science, Google Scholar**112** - 29.
Wang J, Lo K-Y, McCue S, Simpson MJ . 2018 The role of initial geometry in experimental models of wound closing.**Chem. Eng. Sci.**, 221-226. (doi:10.1016/j.ces.2018.01.004) Crossref, Web of Science, Google Scholar**179** - 30.
McCue SW, Jin W, Moroney TJ, Lo K-Y, Chou S-E, Simpson MJ . 2019 Hole-closing model reveals exponents for nonlinear degenerate diffusivity functions in cell biology.**Physica D**, 130. (doi:10.1016/j.physd.2019.06.005) Crossref, Web of Science, Google Scholar**398** - 31.
Callens SJP, Ulyttendaele RJC, Fratila-Apachitei LE, Zadpoor AA . 2020 Substrate curvature as a cue to guide spatiotemporal cell and tissue organization.**Biomaterials**, 119739, 1–22. (doi:10.1016/j.biomaterials.2019.119739) Crossref, PubMed, Web of Science, Google Scholar**232** - 32.
Browning A, Maclaren OJ, Buenzli PR, Lanaro M, Allenby MC, Woodruff MA, Simpson MJ . 2021 Model-based data analysis of tissue growth in thin 3D printed scaffolds.**J. Theor. Biol.**, 110852. (doi:10.1016/j.jtbi.2021.110852) Crossref, PubMed, Web of Science, Google Scholar**528** - 33.
Goriely A . 2017**The mathematics and mechanics of biological growth**. New York, NY: Springer. Crossref, Google Scholar - 34.
Alias A, Buenzli PR . 2017 Modeling the effect of curvature on the collective behavior of cells growing new tissue.**Biophys. J.**, 193-204. (doi:10.1016/j.bpj.2016.11.3203) Crossref, PubMed, Web of Science, Google Scholar**112** - 35.
Alias A, Buenzli PR . 2018 Osteoblasts infill irregular pores under curvature and porosity controls: a hypothesis-testing analysis of cell behaviours.**Biomech. Model Mechanobiol.**, 1357-1371. (doi:10.1007/s10237-018-1031-x) Crossref, PubMed, Web of Science, Google Scholar**17** - 36.
Alias A, Buenzli PR . 2019 A level-set method for the evolution of cells and tissue during curvature-controlled growth.**Int. J. Numer. Meth. Biomed. Engng.**, e3279. (doi:10.1002/cnm.3279) Google Scholar**2019** - 37.
Hegarty-Cremer SGD, Simpson MJ, Andersen TL, Buenzli PR . 2021 Modelling cell guidance and curvature control in evolving biological tissues.**J. Theor. Biol.**, 110658. (doi:10.1016/j.jtbi.2021.110658) Crossref, PubMed, Web of Science, Google Scholar**520** - 38.
Buenzli PR . 2016 Governing equations of tissue modelling and remodelling: a unified generalised description of surface and bulk balance.**PLoS ONE**, e0152582. (doi:10.1371/journal.pone.0152582) Crossref, PubMed, Web of Science, Google Scholar**11** - 39.
Edelstein-Keshet L . 2005**Mathematical models in biology**. Philadelphia, PA: SIAM. Crossref, Google Scholar - 40.
Murray JD . 2002**Mathematical biology: I. An introduction**. New York, NY: Springer. Crossref, Google Scholar - 41.
Li Y, Buenzli PR, Simpson MJ . 2022 Interpreting how nonlinear diffusion affects the fate of bistable populations using a discrete modelling framework.**Proc. R. Soc. A**, 20220013. (doi:10.1098/rspa.2022.0013) Link, Google Scholar**478** - 42.
Shigesada N, Kawasaki K . 1997**Biological invasions: theory and practice**. Oxford, UK: Oxford University Press. Crossref, Google Scholar - 43.
Abramowitz M, Stegun IA (eds). 1972*Handbook of mathematical functions with formulas, graphs, and mathematical tables*, 10th printing, Applied Mathematics Series 55. Washington DC, New York: United States Department of Commerce, National Bureau of Standards. Google Scholar - 44.
Lax PD . 1973**Hyperbolic systems of conservation laws and the mathematical theory of shock waves**. Philadelphia, PA: SIAM. Crossref, Google Scholar - 45.
Taylor CM, Hastings A . 2005 Allee effects in biological invasions.**Ecol. Lett.**, 895-908. (doi:10.1111/j.1461-0248.2005.00787.x) Crossref, Web of Science, Google Scholar**8** - 46.
Maselko J, Showalter K . 1989 Chemical waves on spherical surfaces.**Nature**, 609-611. (doi:10.1038/339609a0) Crossref, Web of Science, Google Scholar**339** - 47.
Cheer A, Vincent JP, Nuccitelli R, Oster G . 1987 Cortical activity in vertebrate eggs I: the activation waves.**J. Theor. Biol.**, 377-404. (doi:10.1016/S0022-5193(87)80217-5) Crossref, PubMed, Web of Science, Google Scholar**124** - 48.
Białecki S, Nałe¸cz-Jawecki P, Kaźmierczak B, Lipniacki T . 2020 Traveling and standing fronts on curved surfaces.**Physica D**, 132215. (doi:10.1016/j.physd.2019.132215) Crossref, Web of Science, Google Scholar**401** - 49.
Westdickenberg MG . 2021 On the metastability of the 1-d Allen–Cahn equation.**J. Dyn. Differ. Equ.**, 1853-1879. (doi:10.1007/s10884-020-09874-z) Crossref, Web of Science, Google Scholar**33** - 50.
Berestycki H, Hamel F . 2002 Front propagation in periodic excitable media.**Comm. Pure App. Math.**, 949-1032. (doi:10.1002/cpa.3022) Crossref, Web of Science, Google Scholar**55** - 51.
El-Hachem M, McCue SW, Jin W, Du Y, Simpson MJ . 2019 Revisiting the Fisher–Kolmogorov–Petrovsky–Piskunov equation to interpret the spreading-extinction dichotomy.**Proc. R. Soc. A**, 20190378. (doi:10.1098/rspa.2019.0378) Link, Google Scholar**475** - 52.
McCue SW, El-Hachem M, Simpson MJ . 2021 Traveling waves, blow-up, and extinction in the Fisher–Stefan model.**Stud. Appl. Math.**, 1-24. (doi:10.1111/sapm.12465) Google Scholar**2021**