Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
Open AccessResearch articles

Curvature dependences of wave propagation in reaction–diffusion models

Pascal R. Buenzli

Pascal R. Buenzli

School of Mathematical Sciences, Queensland University of Technology (QUT), Brisbane, Queensland 4000, Australia

[email protected]

Contribution: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

and
Matthew J. Simpson

Matthew J. Simpson

School of Mathematical Sciences, Queensland University of Technology (QUT), Brisbane, Queensland 4000, Australia

Contribution: Conceptualization, Investigation, Methodology, Writing – review & editing

Google Scholar

Find this author on PubMed

Published:https://doi.org/10.1098/rspa.2022.0582

    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 [110].

    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

    v=cDκ,1.1
    at leading order in the diffusivity D>0 as D0, where κ is the front’s mean curvature [1,2,4].

    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 [1317]. 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,1821]. 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,2326]. 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,2732]. In this context, the eikonal equation represents an understanding of a biological growth law based on migratory and proliferative cellular behaviours [3337].

    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 D1/2 [1,11], so the curvature-dependent term in equation (1.1) is only relevant for curvatures of order O(D1/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

    ut=D2u+F(u),2.1
    describing the spatio-temporal evolution of a quantity u(r,t), where r is a position vector and t is time, with constant diffusivity D>0 and general reaction term F(u). The subscript ut denotes u/t, and 2u is the Laplacian of u. The quantity u may represent the density of entities such as a chemical species, particles or cells forming a biological tissue, which undergo diffusive motion and may be generated and eliminated through the source term F(u). While the exact curvature dependence elucidated below is valid in two- and three-dimensional space, our specific examples and further analyses will focus on two-dimensional domains Ω for simplicity. We supplement equation (2.1) with either inhomogeneous Dirichlet boundary conditions or homogeneous Neumann (no flux) boundary conditions:
    u|Ω=u(Dirichlet) ornu|Ω=0(Neumann) ,
    where n is the inward-pointing unit normal vector to the boundary Ω of the domain. We assume the domain to be initially empty and set u(r,0)=0 in the interior of Ω. With Dirichlet boundary conditions, we set u(r,0)=u at the boundary Ω. With Neumann boundary conditions, we set u(r,0)=uδΩ(r), where δΩ is a surface delta distribution [38]. Formally, u is zero everywhere except at the boundary Ω where u is infinite, such that dru(r,0)=Ωdσ=u|Ω| is proportional to the perimeter of Ω.

    A source term F(u) of particular interest is the logistic growth term F(u)=λu(1u) for a normalized density u[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)=λu and F(u)=λ(1u). 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)=λu(ua)(1u) 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(r,t) is defined to be the time-dependent set of points r in the domain Ω such that u(r,t)=uc for a constant value uc (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 r of the domain has a density u(r,t)=uc for a certain allowable value uc of the solution. The velocity of a contour uc at r can be determined as in level-set methods [13,14]. In short, given a parametrization sr(s,t) of a contour line at time t, the velocity of a point r(s,t) of the contour at fixed s is rt. Differentiating u(r(s,t),t))=uc with respect to t and using the fact that the unit normal to level sets of u in the decreasing direction of u is n=u/|u|, one obtains

    ut+urt=ut|u|nrt=ut|u|v=0,
    where v=nrt is the normal speed in the direction of lower values of u, by definition of n. Given a solution u(r,t), the normal speed of travelling fronts can thus be estimated at any point of the domain and any time by evaluating [5]:
    v=ut|u|=1|u|(D2u+F(u)).2.2
    It is clear from equation (2.2) that this velocity is, in general, space and time dependent and that there may be contributions due to both the diffusion term and the reaction term of equation (2.1), as illustrated in figure 2. Equation (2.2) also describes the speed of travelling-wave solutions that may establish in an infinite domain under appropriate boundary conditions in the limit t. A travelling-wave solution is such that u(r,t)=u0(rct), where c is a constant velocity vector and u0 is the wave profile. For such travelling-wave solutions, ut=u0c and u=u0, so that
    v=ut|u|=u0c|u0|=nc=c,
    where n is the unit vector normal to contours of u0 and c is the normal velocity of the travelling wave. For source terms F(u) and initial and boundary conditions that lead equation (2.1) to establish travelling-wave solutions in the long-time limit, equation (2.2) is such that limtv=nc=c. The value of working with equation (2.2) is that it provides more practical insight that could be relevant in a finite time, finite space experiment where, strictly speaking, travelling-wave solutions never arise.
    Figure 1.

    Figure 1. Numerical simulations of the Fisher–KPP reaction–diffusion equation (F(u)=λu(1u)) on domains of different shapes (white) and linear size L, showing that the speed of travelling fronts u(r,t)=uc (red) is affected by the front’s curvature; fronts uc=0.5 are shown every 1.5 time units, except for the inward-moving fronts of the cross and of the petal that are shown every 0.7 time units. (a) Inward-moving travelling fronts with λ=1, D=0.005; (b) outward-moving travelling fronts with λ=1, D=0.02. In all these simulations, u(r,0)=0 in the domain and u(r,t)=1t on the domain boundary (black). (Online version in colour.)

    Figure 2.

    Figure 2. The curvature of contours u=uc (solid black lines) influences contour propagation speed via both (a) reaction rates and (b) diffusion. (a) The probability of particles (black dots) undergoing a reaction is enhanced at concavities (κ<0) because particles are more likely to encounter there as they accumulate by autocatalytic production, and it is reduced at convexities (κ>0) because particles are less likely to encounter there, as indicated by the blue-shaded areas. (b) Particles (black dots) evenly distributed on a constant density contour u(r,t)=uc at time t move on average with diffusive flux J=Du normal to the contour in the direction of lower density (blue arrows). Near concavities, particles move closer to each other, which increases their density (u>uc). Near convexitites, particles move away from each other, which decreases their density (u<uc). The constant density contour u(r,t+Δt)=uc at time t+Δt is thereby pushed further ahead near concavities than near convexities. (Online version in colour.)

    (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(r,t)=uc 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(r,t)=uc. The normal component of the Laplacian is unn=nTH(u)n=2u/n2, where H(u)=u is the Hessian matrix of u and /n=n. The transverse component of the Laplacian, 2uunn, is always proportional to the mean curvature κ=n of the travelling front (appendix A), i.e. 2uunn=κ|u|, so that

    2u=unnκ|u|.2.3
    Substituting equation (2.3) into equation (2.2) elucidates an exact dependence between the normal velocity of travelling fronts and their mean curvature:
    v(κ)=wDκ,2.4
    where
    w1|u|[Dunn+F(u)]2.5
    is a remaining contribution to v that only depends on the local profile of the solution u(r,t) in the normal direction. Indeed, |u|=u/n and unn are derivatives of u in the direction n.

    We note here that our sign conventions for the unit normal vector n and curvature κ=n are such that κ<0 at concave portions of the domain boundary and κ>0 at convex portions of the domain boundary. This sign convention and interpretation of convexity carries over to travelling fronts u(r,t)=uc inside the domain by interpreting regions where u>uc as occupied, and regions where u<uc 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(κ) 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

    w1D=1|ux|[Duxx+F(u)]tc,2.6
    where convergence to c holds if the solution becomes a travelling wave in the long-time limit, which typically depends on initial and boundary conditions. The term Dκ in equation (1.1) is thus not only valid to O(D2) as D0, but is an exact dependence, which holds for any value of D. While we may have ww1Dc in some asymptotic limits [1,6], we will show below that w in equation (2.4) may also differ significantly from the one-dimensional travelling front speed w1D and travelling-wave speed c due to the proximity of boundaries and other geometric constraints affecting the behaviour of the solution in confined spaces. The contribution w in equation (2.5) only involves the profile of u(r,t) in the normal direction, but this profile is still coupled via the evolution equation with the solution profile in transverse directions. In other words, w may still have implicit dependences upon the curvature κ of the travelling front (which is a transverse property of the front), and these dependences may arise from both reaction and diffusion terms. In the remainder of this article, we investigate these specific contributions in more detail, through numerical simulations, asymptotic behaviours and analytic solutions for a range of domain shapes, boundary conditions and source terms F(u).

    (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 Ω, we discretize a square domain of size L×L large enough to contain Ω using a regular grid with constant space steps Δx=Δy=h=L/(N1). We define

    xi=xmin+ih,i=0,,N1
    and
    yj=ymin+jh,j=0,,N1,
    where xmin=L/2 and ymin=L/2 are offsets allowing us to centre the computation domain around the origin. We discretize time as tk=kΔt, k=0,1,2,, so that tk+1=tk+Δt. Letting uijku(xi,yi,tk) be the numerical approximation of the solution at position (xi,yi) and time tk, the numerical solution of equation (2.1) is stepped in time by solving
    uijk+1=uijk+DΔth2([ui1,jk2uijk+ui+1,jk]+[ui,j1k2uijk+ui,j+1k])+ΔtF(uijk)2.7
    iteratively in k for all indices (i,j) for which (xi,yj)Ω. Dirichlet boundary conditions are implemented by setting
    uijk=u(i,j)such that (xi,yj)Ω,k.
    Neumann boundary conditions are only implemented in this work when Ω is a square domain, or a circular domain. For a square domain, Neumann boundary conditions are implemented using standard second-order discretization across the boundary, so that for points along the boundaries i=0, i=N1, j=0 and j=N1, equation (2.7) is used by defining
    u1,jku1,jk,uN,jkuN2,jk,ui,1kui,1k,ui,Nkui,N2k.
    The initial condition is implemented by setting
    uij0={0,(i,j)such that (xi,yi)Ωu,otherwise.
    For Neumann boundary conditions, this represents the fact that u(r,0)=uδΩ(r) with dru(r,0)=u|Ω|.

    (i) Circular domain

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

    ut=D1r(rur)r+F(u)=D(urr+1rur)+F(u).2.8
    It should be emphasized that this equation is well behaved as r0. Indeed, by symmetry,
    limr0ur=0,2.9
    so that using Bernoulli–L’Hôpital’s rule, we have
    limr0urr=limr0urr.2.10
    To solve equation (2.8) numerically, we discretize the radial axis as ri=ih, i=0,,N1, with h=R/(N1) use second-order finite difference for all spatial derivatives involved and solve
    uik+1=uik+{2DΔth2(ui1k2uik+ui+1k)+ΔtF(uik),i=0,DΔth2(ui1k2uik+ui+1k+12i[ui+1kui1k])+ΔtF(uik),i=1,2,,N1,2.11
    iteratively in k. The time-stepping rule for i=0 is based on equation (2.10), and u1k is set to u1k to satisfy the Neumann-like symmetry condition (2.9) at the origin.

    Dirichlet boundary conditions at the circular boundary Ω are implemented by setting uNk=u. Neumann boundary conditions are implemented by setting uNk=uN2k in equation (2.11).

    (ii) Curvature and velocity

    Numerical estimates of the mean curvature κ are calculated from:

    κ=u|u|=uxxuy22uxuyuxy+uyyux2[ux2+uy2] 3/2,2.12
    using second-order accurate centred finite differences for all first-order and second-order derivatives involved [13]. Numerical estimates of the normal velocity are based on equation (2.2), which involves 2u and |u|=ux2+uy2. These derivatives are also estimated using second-order accurate centred finite differences.

    (iii) Domain boundaries

    Figure 1 illustrates inward-moving fronts and outward-moving fronts on domains Ω 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(θ)=sin4(θ)+cos4(θ);

    (iv)

    A petal with three branches (referred to below as three petal), parametrized in polar coordinates by the radial curve

    R(θ)=R0+acos(nθ),where R0=0.65,a=0.35,n=3.

    Inward-moving front simulations use a computational domain [1,1]×[1,1] (L=2). Outward-moving front simulations use a computational domain [8,8]×[8,8] (L=16). Dirichlet boundary conditions are implemented for all the aforementioned domain shapes. Neumann boundary conditions are only implemented for the circular domain (in polar coordinates) and for the square domain (in Cartesian coordinates). Space and time units are arbitrary.

    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,

    F(u)=λu(1u).3.1
    In the infinite one-dimensional space, this model leads to the establishment of travelling-wave solutions progressing with speed
    c=2Dλ,3.2
    provided the initial profile u(x,0) has compact support [39,40,42].

    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=1. The solution exhibits travelling fronts progressing inward into the circular pore space (figure 3a). Contour lines u=uc are shown every 0.1 increment of uc. The contour line uc=0.5, corresponding to locations with maximum growth rate, is shown in red.

    Figure 3.

    Figure 3. Curvature dependence of front speed in a circular pore (F(u)=λu(1u), λ=1, D=0.005, h=0.01, Δt=0.001 and Dirichlet boundary conditions u|Ω=1). (a) Snapshots of the solution u(r,t) at different times. Contours u=uc are shown every 0.1 increment between uc=0 and uc=1 in dashed white lines, with uc=0.5 shown in solid red. (b) Front speeds v(κ) determined numerically in the whole domain, coloured by the corresponding value of uc (pluses). The asymptotic expression v(κ)cDκ from equation (1.1) is shown as reference (green). (Online version in colour.)

    Numerical estimates of front travelling speed and front curvature obtained at each (discretized) location r of the domain are shown as data points (v(r,t),κ(r,t)) in figure 3b. These data points are coloured by the value uc of the corresponding contour. At all times, there is a well-defined linear correlation v(κ) 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 κ 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 κ 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(κ) 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 κ 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 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 4a. This comparison confirms that the asymptotic growth law (1.1) in figure 3b is well satisfied when and where the solution profile in the normal direction matches the one-dimensional travelling-wave profile. Figure 4b also shows numerical results of inward progressing fronts of the two-dimensional Fisher–KPP model obtained with Neumann boundary conditions at |r|=1.

    Figure 4.

    Figure 4. Comparison between solution profiles u(x,0,t) in a y=0 cross section of the circular domain (inward wave propagation) with either (a) Dirichlet boundary conditions (u|Ω=1) or (b) Neumann boundary conditions. Profiles are shown every 0.5 time units (black). Three phases are distinguished: (1) an establishment phase; (2) a travelling-wave-like phase, where solution profiles are similar to the one-dimensional travelling-wave solution profile (magenta), particularly at low values of u (solid magenta line); (3) a dip-filling phase. The radially symmetric solutions u(x,y,t) are computed numerically in polar coordinates on a domain r[0,1]. The one-dimensional travelling-wave solution profiles (magenta) are a mid-domain snapshot of the solution computed numerically in Cartesian coordinates on the domain x[8,8] with u(8,t)=1, u(8,t)=0, suitably shifted along x. (Online version in colour.)

    (a) Wave propagation phases

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

    (1)

    An establishment phase, where solution profiles are strongly influenced by the proximity of boundaries and by the type of the boundary condition implemented;

    (2)

    A travelling-wave-like phase, where solution profiles normal to the fronts closely match the travelling-wave profiles of the corresponding infinite, one-dimensional space solution;

    (3)

    A dip-filling phase, where fronts coming from opposite ends of the domain and travelling in opposite directions, collide and interact to build up the quantity u symmetrically.

    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 5b, 6b and 7b show that there is a clear correlation between v and κ 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 κ 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(κ) matches the asymptotic relationship (1.1) well, and a dip-filling phase, characterized by a linear relationship v(κ) but with a time-dependent slope.

    Figure 5.

    Figure 5. (a and b) Curvature dependence of front speed in a square pore (F(u)=λu(1u), λ=1, D=0.005, h=0.01, Δt=0.001 and Dirichlet boundary conditions u|Ω=1). Plot specifications are as shown in figure 3. (Online version in colour.)

    Figure 6.

    Figure 6. (a and b) Curvature dependence of front speed in a cross-shaped pore with negative and positive curvature (F(u)=λu(1u), λ=1, D=0.005, h=0.01, Δt=0.001 and Dirichlet boundary conditions u|Ω=1). Plot specifications are as shown in figure 3. (Online version in colour.)

    Figure 7.

    Figure 7. (a and b) Curvature dependence of front speed in a three-petal-shaped pore with negative and positive curvature (F(u)=λu(1u), λ=1, D=0.005, h=0.01, Δt=0.001 and Dirichlet boundary conditions u|Ω=1). Plot specifications are as shown in figure 3. (Online version in colour.)

    Figure 8.

    Figure 8. Solution profiles u(x,0,t) in a y=0 cross section of (a) the circular pore in figure 3; (b) the square pore in figure 5; (c) the cross-shaped pore in figure 6; and (d) the three-petal-shaped pore in figure 7. Profiles are shown every 0.5 time units (black) as well as at the times selected in figures 37 (red). The intersection of the solution profiles with the uc=0.5 level (red dashed line) corresponds to the travelling front shown as the red contour in figures 37. (Online version in colour.)

    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 u0 to u1, 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=uc [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(κ) 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 9b).

    Figure 9.

    Figure 9. Curvature dependence of front speed in a circular pore of radius 4 (F(u)=λu(1u), λ=1, D=0.005, h=0.01, Δt=0.001 and u|Ω=1). (a) Solution profiles u(x,0,t) are shown every 1 time unit (black) and at t=7.0,15.0,20.0 (red). (b) Establishment phase: front speeds differ significantly from the asymptotic expression v(κ)cDκ from equation (1.1) (green). (c and d) Travelling-wave phase: front speeds match the asymptotic expression v(κ)cDκ from equation (1.1) well, except at points far from the wavefront, i.e. near boundaries of the domain where u1 (yellow pluses), and near the centre of the domain where curvature is high. (Online version in colour.)

    The normal velocity v is not always a sole function of curvature. The spread of points (κ,v) in figures 6b and 7b 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 κ 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=uc 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(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. 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 37 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(κ) 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(κ) arise that reinforce this tendency.

    To analyse these additional contributions to v(κ) in the dip-filling phase, we now assume that the solution is approximately radially symmetric near the centre of the domain r0. Since n points in the opposite direction to the radial axis, /n=n=/r. The contribution to v(κ) due to diffusion in the normal direction captured by the term Dunn/|u| in equation (2.5) becomes, using equations (2.9) and Bernoulli–L’Hôpital’s rule (2.10):

    Dunn|u|DurrurDur/rurDr,r0.
    Since κ=1/r, we obtain Dunn/|u|Dκ as r0, so that the speed of travelling fronts in equation (2.4) becomes
    v(κ)2DκF(u)un,r0 or κ.3.3
    We conclude that in the dip-filling phase, both the normal component of diffusion and the transverse component of diffusion contribute a term Dκ each to front speed. While this is valid at all times for radially symmetric solutions as κ, it is only valid as t for non-circular domains since moving fronts in such domains evolve into circles only in the long-time limit. It is interesting to note here that some analyses of flame propagation also suggest flame displacement speed to be proportional to curvature with a factor twice the thermal diffusivity [8]. In figure 10, we illustrate this extra contribution of curvature to front speed in the dip-filling phase by presenting numerical simulations of the diffusion equation (F(u)=0) in a circular domain with Dirichlet boundary conditions u=1. These simulations show that without reaction terms, front speed in the dip-filling phase indeed converges to the time-independent linear relationship v(κ)2Dκ. The speed of a travelling front, as it progresses in space, is time dependent due its changing curvature. The diffusion equation is an example for which the solution does not develop into a travelling wave in infinite space.
    Figure 10.

    Figure 10. Curvature dependence of front speed in a circular pore for the diffusion equation with Dirichlet boundary conditions (F(u)=0, D=0.005, h=0.01, Δt=0.001 and u|Ω=1). (a) Solution profiles are shown every 2.5 time units (black) as well as at t=3.5, t=10.0 and t=100.0 (red). (b) At early time (t=3.5) and small curvature (far from the centre), the relationship v(κ) has a slope close to D, which transitions into the asymptotic relationship v(κ)2Dκ at large curvature (see equation (3.3)). (c and d) With time, all fronts of the solution move with the time-independent speed v(κ)=2Dκ. (Online version in colour.)

    Importantly, the contribution to v(κ) 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 ur/rurr as r0 and un=ur, we have

    1un(r,t)1rurr(0,t)=κunn(0,t),r0 or κ.3.4
    Thus, in the dip-filling phase:
    v(κ)2DκF(u(0,t))unn(0,t)κ,κ.3.5
    Equation (3.5) shows that in the dip-filling phase, fronts evolve with a normal velocity dominantly proportional to curvature. The proportionality coefficient has an explicit, constant contribution due to diffusion, and a time-dependent contribution related to the autocatalytic evolution of the solution at the centre of the dip (r=0) governed by the reaction term F. It should be emphasized that the second term on the right-hand side of equation (3.5) may still contain implicit dependences upon D via u(0,t) and unn(0,t).

    (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(κ) 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(κ); and (ii) how far from the origin r0 we can expect the asymptotic behaviour of v(κ) 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 D0 in the dip-filling phase. We show that in this regime, the eikonal equation v(κ) 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(κ) observed at long times in figures 3 and 57, and captured in equation (3.5). First, let us show that the velocity field v(r,t) of all moving fronts is independent of time if D=0. From equation (2.2) with D=0, we have

    vt(r,t)=F(u)ut|u|F(u)|u|2(|u|)t=F(u)F(u)|u|F(u)|u|2F(u)|u|=0,3.6
    where we have used
    |u|t=tuu=utu|u|=F(u)uu|u|=F(u)|u|.
    In the dip-filling phase where travelling fronts become circular, the curvature of a front of radius r is 1/r. The curvature of this front increases with time (in absolute value) as the inward-moving front comes closer to the origin r=0. At a fixed location in space, however, the curvature map κ(r,t) for circular fronts about the origin is simply given by κ(r,t)=1/r, where r=|r|, and it is therefore independent of time. Since both the velocity map v(r,t) and curvature map κ(r,t) are time independent in the zero diffusion limit of the dip-filling phase, the relationship v(κ) is also time independent. We can conclude that the steepening with time of the linear relationship v(κ) in our numerical simulations in figures 37 is entirely due to diffusive effects.

    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 11a). We observe that the relationship v(κ) 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(κ) no longer evolves in time and is approximately given by v(κ)=γκ, where γ>0 is a constant (figure 11b).

    Figure 11.

    Figure 11. Curvature dependence of front speed when diffusion is turned off at t=9.0 in the two-dimensional analogue of the Fisher–KPP model with Dirichlet boundary conditions (F(u)=λu(1u), λ=1, D=0.005, h=0.01, Δt=0.001 and u|Ω=1). (a) Solution profiles are shown every 0.5 time units (black) and at t=9.0 (red). The parameter γ in equation (3.9) that best matches the solution profile at t=9.0 near the domain centre (solid magenta curve) predicts the slope of the relationship v(κ)=γκ very well (magenta line in (b)). Dashed magenta lines show U(r,t) given by equation (3.9) every 0.5 time units from t=9.5 to t=12.0. The constant A00.838 in equation (3.9) is taken from the numerical solution u(0,9.0). For increased accuracy in (b), curvature κ at the radial coordinate r is taken to be 1/r rather than being estimated numerically from the solution u. (Online version in colour.)

    The proportionality constant γ 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 t0 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>t0 is given by a function U(r,t), where r0 is the radial coordinate. From equation (2.2) and assuming v(κ)=γκ=γ/r, we have:

    v=F(U)Ur=γr,
    or, equivalently,
    UrF(U)=rγ.3.7
    Integrating the differential equation in r (3.7) for U(r,t) can be done explicitly for several choices of the reaction term F(u). For logistic growth F(u)=λu(1u), 1/F(u)=λ1(d/du)[ln(u)ln(1u)], so that
    U(r,t)=11+A(t)exp{λr2/(2γ)},3.8
    where the integration constant A(t)=1/U(0,t)1 contains the only time dependence. Differentiating this solution with respect to t gives
    Ut=AtAU(1U),
    which shows that A(t)=A0exp{λ(tt0)}. Therefore, for t>t0,
    U(r,t)=11+A0exp{λr2/(2γ)λ(tt0)},3.9
    which means that the solution converges exponentially fast to the uniform asymptotic solution limtU(r,t)=1, with contour lines that move inward with velocity v(r)=γ/r.

    The value of γ that best fits the numerical solution profile at t=9.0 around r=0, as shown in figure 11a, is γ0.047. This value of γ found from the solution profile predicts the slope of the growth law v(κ) very well, see figure 11b. The solution in equation (3.9) differs from the numerical solution at locations where fronts have small curvature (|x|0.2, corresponding to |κ|5, see figure 11a). At these curvatures, numerical estimates of v(κ) in figure 11b deviate from the linear relationship v(κ)=γκ. Figure 11a 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(κ) 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(κ) by integrating, from equation (2.4),

    UrF(U)=1v(1/r)2D/r,
    with respect to r, e.g. using numerical quadrature. This correspondence may be useful to predict the speed of travelling fronts in the dip-filling phase from a single snapshot in time of the solution profile, such as that provided by a biological experiment. On the other hand, measuring curvature-dependent speeds of travelling fronts may allow the estimation of solution profiles that may otherwise be difficult to measure, such as in wildfires.

    In summary, the dip-filling phase is characterized by travelling front velocities v(κ) that are affected by curvature due to several contributions: (i) the transverse component of diffusion, always exactly contributing a term Dκ; (ii) the normal component of diffusion, contributing an extra term Dκ at large curvature κ; and (iii) reaction terms, which provide additional contributions to the normal speed that also becomes linear in κ at large curvature κ, 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(κ), 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)=λu, and a saturated growth model, where F(u)=λ(1u). 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

    F(u)=λu.3.10
    The solution can be found by using the Fourier–Laplace transform of u(r,t) in time and solving the resulting eigenvalue problem of the Laplacian in polar coordinates with radial symmetry. The general solution is a superposition of Bessel functions of the first and second kind [43]. Only Bessel functions of the first kind, Jn(x), are well defined at x=0, and only the angular mode n=0 is radially symmetric, so solutions have the form
    u(r,t)=[a0+s=1aksJ0(ksr)exp{ks2Dt}]exp{λt},
    where a0 and aks, s=1,2,, are constants that depend on the initial profile u(r,0). The radial modes ks are determined by the boundary condition at r=R. If λ=0, the expression mentioned earlier is the general solution of the diffusion equation with radial symmetry. In this case, we can satisfy the Dirichlet boundary condition u(R,t)=u=1 by setting a0=1 and ks such that ksR=j0,s, where j0,s>0, s=1,2,, are the roots of the Bessel function J0. If λ>0, it is more consistent to satisfy a time-dependent Dirichlet boundary condition u(R,t)=exp{λt}, in which case a0=1 and ks are as mentioned earlier. Alternatively, given that J0(x)=J1(x), we can enforce Neumann boundary conditions by selecting modes ks such that ksR=j1,s, where j1,s>0, s=1,2,, are the roots of the Bessel function J1.

    To solve the equivalent problem with the saturated growth reaction term

    F(u)=λ(1u),
    we let w=1u, so that wt=D2wλw is Skellam’s model with the opposite sign for λ. The solution u=1w with Dirichlet boundary conditions is thus given by
    u(r,t)=1s=1aksJ0(ksr)exp{ks2Dt}exp{λt}.
    At leading order in the long-time limit t, we thus have, for the diffusion model (udiff), Skellam’s model with Dirichlet (uSkellamD) and Neumann (uSkellamN) boundary conditions, and for the saturated growth model (usat):
    udiff(r,t)1ak1J0(k1r)exp{k12Dt},(F(u)=0,u(1,t)=1),3.11
    uSkellamD(r,t)[1ak1J0(k1r)exp{k12Dt}]exp{λt},(F(u)=λu,u(1,t)=exp{λt}),3.12
    uSkellamN(r,t)[a0ak1J0(k1r)exp{k12Dt}]exp{λt},(F(u)=λu,ur(1,t)=0)3.13
    andusat(r,t)1ak1J0(k1r)exp{k12Dt}exp{λt},(F(u)=λ(1u),u(1,t)=1),3.14
    where the signs in front of a0,ak1 and ak1 are set such that these constants are all positive, and where we assume R=1, so that k1=j0,12.405 and k1=j1,13.832. We can use the explicit long-term solutions in equations (3.11)–(3.14) to estimate the speed of moving fronts in the dip-filling phase according to equations (2.4)–(2.5). By using /n=/r, we have
    Dunnun=Durrurk1DJ0(k1r)J0(k1r),t,
    (with k1 replaced by k1 for Skellam’s model with Neumann boundary conditions). From Bessel’s equation satisfied by J0, J0(x)/J0(x)+1/x+J0(x)/J0(x)=0, so that by substituting x=k1r and using κ=1/r and J0(x)=J1(x) [43], we obtain
    DunnunDκ+k1DJ0(k1/κ)J1(k1/κ),t,
    in all four models. The velocity of fronts is thus given by
    v(κ)k1DJ0(k1/κ)J1(k1/κ)+F(u)ur,t,2Dκ+F(u)ur,t,κ,3.15
    where the last asymptotic expression uses J0(x)=1x2/4+O(x4), and J1(x)=x/2+O(x3) as x0 [43]. This result is an explicit verification of the asymptotic behaviour derived in equation (3.3). Equation (3.15) also provides the contribution to v(κ) due to diffusion for any curvature in the long-time limit, but the difference with the large curvature expression is very small (figure 10).

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

    vSkellamD(κ)(k1Dλk1)J0(k1/κ)J1(k1/κ)+λexp{k12Dt}k1ak1J1(k1/κ),t,2Dκ2λk12(exp{k12Dt}ak11)κ,t,κ.

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

    vsat(κ)(k1D+λk1[1+J(κ)e(k22k12)Dt])J0(k1/κ)J1(k1/κ),t,2Dκ2λk12(1+ak2ak1(1k22k12)e(k22k12)Dt)κ,t,κ,
    where
    J(κ)=ak2ak1(J0(k2/κ)J0(k1/κ)k2J1(k2/κ)k1J1(k1/κ))ak2ak1(1k22k12),κ.
    It can also be checked explicitly from equations (3.11)–(3.14) that the reaction-dependent expressions obtained earlier by calculating F(u)/ur in the limit r0 correspond exactly to [F(u(0,t))/urr(0,t)]κ, as predicted by equation (3.5).

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

    limtvSat(κ)=2Dκ2λk12κ,3.16
    as confirmed by our numerical simulations (figures 12 and 13). Figure 13 shows that the front speed law v(κ) in the two-dimensional Fisher–KPP model also converges to equation (3.16) in the long-time limit, which is consistent with the fact that the saturated growth model is the linearization of the Fisher–KPP model about the steady state u=1. At early times, however, the Fisher–KPP model is similar to Skellam’s model. The slope of v(κ) therefore first increases in time before converging to that of equation (3.16) in all pore shapes (figures 3 and 57).
    Figure 12.

    Figure 12. Curvature dependence of front speed in a circular pore for the saturated growth model with Dirichlet boundary conditions (F(u)=λ(1u), λ=1, D=0.05, h=0.01, Δt=0.0001 and u|Ω=1). (a) Solution profiles are shown every 0.15 time units (black) as well as at the times t=0.06,2.25,4.50 (red). (b)–(d) With time, all fronts of the solution move with a speed approaching the asymptotic limit vSat(κ) in equation (3.16) (blue). In these simulations, diffusivity is increased by a factor 10 compared with previous inward-moving front simulations, to prevent the reaction term from dominating the evolution in the dip-filling phase. (Online version in colour.)

    Figure 13.

    Figure 13. Slope of the linear relationship v(κ) in the dip-filling phase versus time. The two-dimensional analogue of the Fisher–KPP model, the saturated growth model, and the strong Allee model all converge to the linear relationship v(κ) given by equation (3.16). In these simulations, all models are solved on the circular domain with Dirichlet boundary conditions and D=0.05, h=0.05 and Δt=0.001. (Online version in colour.)

    (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)=λu(ua)(1u) with 0<a<1. This model supports travelling waves with speed

    c=2Dλ(12a),
    in one-dimensional infinite space under appropriate boundary and initial conditions [4]. In higher dimensions, the invasion of available space by the solution depends on the size and the shape of the initial condition [4,41]. Numerical simulations of this model are shown in figure 14 with a=2/5 and λ=1/(1a)=5/3 so that F(u)1u as u1 and c>0. Here too, the relationship v(κ) converges to that of the saturated growth model in equation (3.16) as t (figure 13).
    Figure 14.

    Figure 14. Curvature dependence of front speed in a circular pore for the strong Allee growth model with Dirichlet boundary conditions (F(u)=λu(ua)(1u) with a=2/5, λ=1/(1a)=3/5, D=0.05, h=0.01, Δt=0.0001 and u|Ω=1). For this model, the minimum speed of travelling waves in infinite one-dimensional space is c=2Dλ(1/2a). (a) Solution profiles are shown every 0.15 time units (black) as well as at the times t=0.45,2.25and9.0 (red). (b)–(d) With time, all fronts of the solution move with a speed approaching the asymptotic limit vSat(κ) in equation (3.16) (magenta). The asymptotic expression v(κ)cDκ from equation (1.1) is also shown for reference (green). Diffusivity is the same as in the saturated growth model for comparison in figure 13. (Online version in colour.)

    Our numerical observation that the linear relationship v(κ) 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(κ) depends on urr(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κ 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=cDκ 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κ 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)=λu), the saturated growth model (F(u)=λ(1u)), as well as the strong Allee source term (F(u)=λu(ua)(1u)), 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 κ that arises from transverse diffusion, and the fact that mean curvature flows evolve interfaces into shrinking circles [1821]. 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(1u2) [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 κ=n. The mean curvature here is defined as the sum of principal curvatures, rather than their average [13,24], and n=u/|u| is the unit normal vector to level sets of u pointing towards decreasing values of u. We start by expanding the divergence of n using standard rules of differentiation:

    κ=n=(u|u|)=2u|u|u(1|u|)=2u|u|+u(|u|)|u|2.A 1
    Furthermore,
    (|u|)=uu=(u)Tu+uTu2|u|=H(u)u|u|,A 2
    where
    H(u)=u={2uxixj}i,j=1m
    is the Hessian matrix collecting all second partials of u in m dimensions of space (m=2,3). We thus have, from equations (A 1) and (A 2):
    κ=2u|u|+1|u|3uH(u)u=1|u|(2unTH(u)n).A 3
    By defining the derivative in the normal direction /n as /n=n, we also see that
    unn=2un2=n(nu)=n(|u|)=nTH(u)n,
    where we used nu=|u| and equation (A 2). Reorganizing equation (A 3), the Laplacian can therefore be decomposed into a normal component unn=nTH(u)n, and a transverse component proportional to the mean curvature:
    2u=unnκ|u|.A 4

    Footnotes

    Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

    References