Curvature dependences of wave propagation in reaction-diffusion models

Reaction-diffusion waves in multiple spatial dimensions advance at a rate that strongly depends on the curvature of the wave fronts. 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 wave front. Inward-propagating waves are characterised 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.


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][2][3][4][5][6][7][8][9][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 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 of travelling fronts is given by at leading order in the diffusivity > 0 as → 0, where is the front's mean curvature [1,2,4]. In principle, Eq. (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 reactiondiffusion 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][14][15][16][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 Eq. (1), round off the initial front shape in such a way that the front becomes circular before shrinking into a point [13,[18][19][20][21]. Such rounding of initial shape is commonly observed in experimental systems [3,7,9,22]. Many variants of mean curvature flows, such as lengh-preserving flows, area-preserving flows, and Willmore flows are proposed as effective models of the evolution of fronts or interfaces [13,14,21,[23][24][25][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][28][29][30][31][32]. In this context, the eikonal equation represents an understanding of a biological growth law based on migratory and proliferative cellular behaviours [33][34][35][36][37].
For the eikonal equation (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) 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 is proportional to 1/2 [1,11], so the curvature-dependent term in Eq. (1) is only relevant for curvatures of order O( −1/2 ). This corresponds to slowly moving highly curved fronts, such as where circular wave fronts 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 Eq. (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 . 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 Eq. (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 Eq. (1) is an exact, non perturbative contribution that originates from the transverse component of diffusion (diffusion parallel to the wave front). Further dependences upon 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.    (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 ( , ) = at time move on average with diffusive flux = − ∇ normal to the contour in the direction of lower density (blue arrows). Near concavities, particles move closer to each other, which increases their density ( > ). Near convexitites, particles move away from each other, which decreases their density ( < ). The constant density contour ( , + Δ ) = at time + Δ is thereby pushed further ahead near concavities than near convexities.

Model description and theoretical developments
We consider the single-species reaction-diffusion model describing the spatio-temporal evolution of a quantity ( , ), where is a position vector and is time, with constant diffusivity > 0 and general reaction term ( ). The subscript denotes / , and ∇ 2 is the Laplacian of . The quantity may represent the density of entities such as a chemical species, particles, or cells forming a biological tissue, that undergo diffusive motion, and that may be generated and eliminated through the source term ( ). 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 Eq. (2) with either inhomogeneous Dirichlet boundary conditions, or homogeneous Neumann (no flux) boundary conditions: where is the inward-pointing unit normal vector to the boundary of the domain. We assume the domain to be initially empty and set ( , 0) = 0 in the interior of . With Dirichlet boundary conditions, we set ( , 0) = * at the boundary . With Neumann boundary conditions, we set ( , 0) = * δ ( ), where δ is a surface delta distribution [38]. Formally, is zero everywhere except at the boundary where is infinite, such that ∫ d ( , 0) = ∫ d = * | | is proportional to the perimeter of . A source term ( ) of particular interest is the logistic growth term ( ) = (1 − ) for a normalised density ∈ [0, 1], in which case Eq. (2) is the two-dimensional analogue of the Fisher-KPP model [11,39]. In this model, = 0 represents the unstable, excitable rest state, and = 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 ( ) leads to travelling-wave solutions, nor that it possesses rest states and excited states. We will state results valid for general source terms ( ), and will illustrate properties of fronts of the solution to Eq. (2) with linear source terms, ( ) = and ( ) = (1 − ). These linear models do not support travelling-wave solutions. The time-dependent speed of their travelling fronts enable us to illustrate important properties of the eikonal equation. We also consider a bistable source term ( ) = ( − )(1 − ) with 0 < < 1, a common model for population growth subject to a strong Allee effect. The rest state = 0 and excited state = 1 in this model are both stable, yet this model still exhibits travelling-wave transitions from = 0 to = 1 or from = 1 to = 0 depending on and on the shape and size of initial conditions [4,41].
Speed of travelling fronts. A travelling front of ( , ) is defined to be the time-dependent set of points in the domain such that ( , ) = for a constant value ( Figure 1). A convenient way to estimate the velocity of travelling fronts is to calculate the velocity of contours of at each location of the domain. Each location of the domain has a density ( , ) = for a certain allowable value of the solution. The velocity of a contour at can be determined as in level-set methods [13,14]. In short, given a parametrisation ↦ → ( , ) of a contour line at time , the velocity of a point ( , ) of the contour at fixed is . Differentiating ( ( , ), )) = with respect to , and using the fact that the unit normal to level sets of in the decreasing direction of is = −∇ /|∇ |, one obtains where = · is the normal speed in the direction of lower values of , by definition of . Given a solution ( , ), the normal speed of travelling fronts can thus be estimated at any point of the domain and any time by evaluating [5]: It is clear from Eq. (3) that this velocity is, in general, space and time dependent, and that there may be contributions due to both the diffusion term and reaction term of Eq. (2), as illustrated in where is the unit vector normal to contours of 0 , and is the normal velocity of the travelling wave. For source terms ( ) and initial and boundary conditions that lead Eq. (2) to establish travelling-wave solutions in the long-time limit, Eq. (3) is such that lim →∞ = · = . The value of working with Equation (3) 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.
Exact curvature dependence due to transverse diffusion. Equation Substituting Eq. (4) into Eq. (3) elucidates an exact dependence between the normal velocity of travelling fronts, and their mean curvature: where is a remaining contribution to that only depends on the local profile of the solution ( , ) in the normal direction. Indeed, −|∇ | = / and are derivatives of in the direction .
We note here that our sign conventions for the unit normal vector and curvature = ∇ · 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 ( , ) = inside the domain by interpreting regions where > as occupied, and regions where < as empty. With this convention, inward-moving circular fronts of the models ( ) we consider have negative curvature, and outward-moving circular fronts have positive curvature.
The expression ( ) in Eq. (5) matches the asymptotic expression (1) from singular perturbation theory if is the speed of the travelling wave in the corresponding one-dimensional problem on an infinite domain, when such waves exist [1,4,6]. The relationship between and can be assessed by comparing Eq. (6) with the speed of travelling fronts in 1D, given from Eq. (3) by where convergence to holds if the solution becomes a travelling wave in the long-time limit, which typically depends on initial and boundary conditions. The term − in Eq. (1) is thus not only valid to O( 2 ) as → 0, but is an exact dependence, which holds for any value of . While we may have ∼ 1D ∼ in some asymptotic limits [1,6], we will show below that in Eq. (5) may also differ significantly from the one-dimensional travelling front speed 1D and travelling-wave speed , due to the proximity of boundaries and other geometric constraints affecting the behaviour of the solution in confined spaces. The contribution in Eq. (6) only involves the profile of ( , ) in the normal direction, but this profile is still coupled via the evolution equation with the solution profile in transverse directions. In other words, 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 the paper, 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 ( ).

Numerical simulations
Numerical simulations of Eq. (2) 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 scheme are possible (such as total variation diminishing Runge-Kutta methods) but curvature effects are known to be more sensitive to spatial discretisation than temporal discretisation [14,36].
For general domain shapes , we discretise a square domain of size × large enough to contain using a regular grid with constant space steps Δ = Δ = ℎ. We define = 0, . . . , − 1, where min = − /2 and min = − /2 are offsets allowing us to centre the computation domain around the origin. We discretise time as = Δ , = 0, 1, 2, . . ., so that +1 = + Δ . Letting ≈ ( , , ) be the numerical approximation of the solution at position ( , ) and time , the numerical solution of Eq. (2) is stepped in time by solving iteratively in for all indices ( , ) for which ( , ) ∈ . Dirichlet boundary conditions are implemented by setting 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 discretisation across the boundary, so that for points along the boundaries The initial condition is implemented by setting otherwise.
For Neumann boundary conditions, this represents the fact that Circular domain. For circular domains with radially symmetric solutions ( , ), where = | |, we implement a FTCS scheme by first expressing Eq. (2) in polar coordinates and assuming circular symmetry. The reaction-diffusion equation becomes It should be emphasised that this equation is well behaved as → 0. Indeed, by symmetry, so that using Bernoulli-L'Hôpital's rule, we have To solve Eq. (9) numerically, we discretise the radial axis as = ℎ, = 0, . . . − 1, use second order finite difference for all spatial derivatives involved, and solve iteratively in . The time-stepping rule for = 0 is based on Eq. (11), and −1 is set to 1 to satisfy the Neumann-like symmetry condition (10) at the origin. Dirichlet boundary conditions at the circular boundary are implemented by setting = * . Neumann boundary conditions are implemented by setting = −2 in Eq. (12).

Curvature and velocity.
Numerical estimates of the mean curvature are calculated from: 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 Eq. (3), which involves ∇ 2 and |∇ | = √︃ 2 + 2 . These derivatives are also estimated using second-orderaccurate centred finite differences. Figure 1 illustrates inward-moving fronts and outward-moving fronts on domains specified by the following boundaries:

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 above. In this model, the reaction term represents logistic growth, In the infinite one-dimensional space, this model leads to the establishment of travelling wave solutions progressing with speed provided the initial profile ( , 0) has compact support [39,40,42].  Numerical estimates of front travelling speed and front curvature obtained at each (discretised) location of the domain are shown as data points ( , ), ( , ) in Figure 3b. These data points are coloured by the value of the corresponding contour. At all times, there is a well defined linear correlation ( ) between front speed and curvature. At = 3.0, several data points with low density match the asymptotic relationship (1), shown in green, reasonably well. Numerical estimates of and or at very low density at = 3.0 may be less accurate as they correspond to points near the centre of the domain, where spatial variations in are small. A value of less than −30 correspond to radii of curvature less than about 3 grid steps ℎ. High density data points are close to the boundary and they deviate significantly from the relationship (1). At times = 7.0 and = 15.0, the relationship ( ) remains mostly linear, but the slope of the relationship is time-dependent, in contrast to Eq. (1).
The results shown in Figure 3 suggest that the asymptotic growth law (1) where the linear dependence of upon has slope − , and where in Eq. (5) is given by the one-dimensional travelling wave speed , Eq. (15), 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) at the same time, or at all. Equations (6)- (7) show that is well approximated by if the solution profile in the normal direction 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 onedimensional travelling-wave profile in Figure 4a. This comparison confirms that the asymptotic growth law (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 condition at | | = 1.  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 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 symmetrically. Numerical simulations of inward progressing fronts of the Fisher-KPP model in the square pore shape, cross pore shape, and 3-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 toward 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 and that becomes increasingly linear as time progresses. However, as before, this relationship is strongly time dependent. The asymptotic relationship (1) between 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 ( ) matches the asymptotic relationship (1) well, and a dip-filling phase, characterised by a linear relationship ( ) but with a time-dependent slope.
The establishment phase results from complex interactions between 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 growth of the density 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 wave front, i.e., the characteristic length over which the solution profile in the normal direction transitions from ≈ 0 to ≈ 1, and the one-dimensional travelling wave speed . The establishment phase transitions into the travelling-wave phase when the solution profile approximates the shape of the one-dimensional travelling wave front 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 4 , the inverse of the steepness of the solution at = [40]. This width is directly related to the travelling wave speed . By decreasing diffusivity, or equivalently, by increasing domain size [9], the travelling-wavelike phase lasts longer, and the match between ( ) and the asymptotic relationship (1) improves ( Figure 9). In 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 wave front 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).
The normal velocity is not always a sole function of curvature. The spread of points ( , ) in Figures 6b, 7b indicates that other detail of the solution may influence the normal velocity. The curves formed by the spread of points in these figures are due to only calculating and at discretised 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 = splits into multiple disconnected fronts ( = 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 ( , ) 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) is well satisfied sufficiently far from the boundaries, but the correction term due to curvature is small and tends to zero as → ∞. 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. 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 ( ) of front speed upon curvature, although with a different slope compared to the asymptotic perturbation result (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 due initially to the explicit curvature-dependence of normal velocity in Eq. (5) [18,19]. As time progresses, further curvature dependences of ( ) arise that reinforce this tendency.

Normal diffusion in the dip-filling phase.
To analyse these additional contributions to ( ) in the dip-filling phase, we now assume that the solution is approximately radially symmetric near the centre of the domain → 0. Since points in the opposite direction to the radial axis, / = · ∇ = − / . The contribution to ( ) due to diffusion in the normal direction captured by the term /|∇ | in Eq. (6) becomes, using Eqs (10), and Bernoulli-L'Hôpital's rule (11): Since = −1/ , we obtain /|∇ | ∼ − as → 0, so that the speed of travelling fronts in Eq. (5) becomes We conclude that in the dip-filling phase, both the normal component of diffusion and the transverse component of diffusion contribute a term − each to front speed. While this is valid at all times for radially symmetric solutions as → −∞, it is only valid as → ∞ 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 ( ( ) = 0) in a circular domain with Dirichlet boundary conditions * = 1. These simulations show that without reaction terms, front speed in the dip-filling phase indeed converges to the time-independent linear relationship ( ) ∼ −2 . 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.
Importantly, the contribution to ( ) 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 perturbations 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.
Reaction-rate dependence in the dip-filling phase. The contribution due to the reaction term in Eq. (16) also possesses a dominant linear dependence upon curvature. Indeed, proceeding similarly to above and using / ∼ as → 0 and = − , we have Thus, in the dip-filling phase: Equation (18) 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 ( = 0) governed by the reaction term . It should be emphasised that the second term on the right hand side of Eq. (18) may still contain implicit dependences upon via (0, ) and (0, ). (16) and (18) show that the dip-filling phase can be highly influenced by reaction terms since fronts of the solution may progress in space due to increasing locally ( Figure 2). However, the reaction-rate-dependent contribution to ( ) on the right hand side of Eq. (18) is still coupled to diffusive effects. In this paragraph, we therefore investigate two remaining important points: (i) the influence of diffusion in the reaction-rate-dependent contribution to ( ); and (ii) how far from the origin → 0 we can expect the asymptotic behaviour of ( ) in Eq. (18) to hold in practice.

Zero-diffusion limit. Equations
To get insights into the velocity of fronts due to reaction terms only, we consider the zero diffusion limit → 0 in the dip-filling phase. We show that in this regime, the eikonal equation ( ) becomes time independent for any choice of ( ). In other words, diffusive effects are entirely responsible for the time dependence of the slope of ( ) observed at long times in where we have used In the dip-filling phase where travelling fronts become circular, the curvature of a front of radius is −1/ . The curvature of this front increases with time (in absolute value) as the inward-moving front comes closer to the origin = 0. At a fixed location in space, however, the curvature map ( , ) for circular fronts about the origin is simply given by ( , ) = −1/ , where = | |, and it is therefore independent of time. Since both the velocity map ( , ) and curvature map ( , ) are time independent in the zero-diffusion limit of the dip-filling phase, the relationship ( ) is also time independent. We can conclude that the steepening with time of the linear relationship ( ) in our numerical simulations in Figures 3-7 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 = 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 ( ) in the dip-filling phase is linear with a slope that becomes steeper with time until = 9.0. After this point, the relationship ( ) no longer evolves in time and is approximately given by ( ) = − , where > 0 is a constant (Figure 11b).
The proportionality constant is directly related to the profile of the solution at = 9.0 and its subsequent evolution at times > 9.0, as follows. Let 0 be the time at which is set to zero. When = 0 and with radial symmetry (dip-filling phase), the solution profile at time > 0 is given by a function ( , ), where 0 is the radial coordinate. From Eq. (3) and assuming ( ) = − = / , we have: or, equivalently, Integrating the differential equation in (20) where the integration constant ( ) = 1/ (0, ) − 1 contains the only time dependence. Differentiating this solution with respect to gives which shows that ( ) = 0 exp{− ( − 0 )}. Therefore, for > 0 , which means that the solution converges exponentially fast to the uniform asymptotic solution lim →∞ ( , ) = 1, with contour lines that move inward with velocity ( ) = / . The value of that best fits the numerical solution profile at = 9.0 around = 0 in Figure 11a, is ≈ 0.047. This value of found from the solution profile predicts the slope of the growth law ( ) very well, see Figure 11b. The solution in Eq. (22) differs from the numerical solution at locations where fronts have small curvature (| | 0.2, corresponding to 5, see Figure 11a). At these curvatures, numerical estimates of ( ) in Figure (11 with respect to , for example 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 characterised by travelling front velocities ( ) that are affected by curvature due to several contributions: (i) the transverse component of diffusion, always exactly contributing a term − ; (ii) the normal component of diffusion, contributing an extra term − at large curvature → −∞; (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.
Linearised models. To understand in more detail how solution profiles converge in the dipfilling phase to profiles that precisely give rise to linear contributions in the growth law ( ), 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 ( ) = , and a saturated-growth model, where ( ) = (1 − ). Skellam's model corresponds to linearising the logistic reaction rate about the unstable rest state = 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 linearising the logistic reaction rate about the stable excited state = 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 Eq. (2) for a radially symmetric solution The solution can be found by using the Fourier-Laplace transform of ( , ) 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, ( ), are well defined at = 0, and only the angular mode = 0 is radially symmetric, so solutions have the form where 0 and , = 1, 2, . . ., are constants that depend on the initial profile ( , 0). The radial modes are determined by the boundary condition at = . If = 0, the expression above is the general solution of the diffusion equation with radial symmetry. In this case, we can satisfy the Dirichlet boundary condition ( , ) = * = 1 by setting 0 = 1 and such that = 0, , where 0, > 0, = 1, 2, . . ., are the roots of the Bessel function 0 . If > 0, it is more consistent to satisfy a time-dependent Dirichlet boundary condition ( , ) = exp{ }, in which case 0 = 1 and are as above. Alternatively, given that 0 ( ) = − 1 ( ), we can enforce Neumann boundary conditions by selecting modes such that = 1, , where 1, > 0, = 1, 2, . . ., are the roots of the Bessel function 1 .
To solve the equivalent problem with saturated growth reaction term is Skellam's model with the opposite sign for . The solution = 1 − with Dirichlet boundary conditions is thus given by At leading order in the long-time limit → ∞, we thus have, for the diffusion model ( diff ), Skellam's model with Dirichlet ( D Skellam ) and Neumann ( N Skellam ) boundary conditions, and for the saturated growth model ( sat ): where the signs in front of 0 , 1 and 1 are set such that these constants are all positive, and where we assume = 1, so that 1 = 0,1 ≈ 2.405, and 1 = 1,1 ≈ 3.832. We can use the explicit long-term solutions in Eq. (24)- (27) to estimate the speed of moving front in the dip-filling phase according to Eq. (5)- (6). Using / = − / , we have , → ∞, (with 1 replaced by 1 for Skellam's model with Neumann boundary conditions). From Bessel's equation satisfied by 0 , 0 ( )/ 0 ( ) + 1/ + 0 ( )/ 0 ( ) = 0, so that substituting = 1 and using = −1/ and 0 ( ) = − 1 ( ) [43], we obtain in all four models. The velocity of fronts is thus given by where the last asymptotic expression uses 0 ( ) = 1 − 2 /4 + O( 4 ), and 1 ( ) = /2 + O( 3 ) as → 0 [43]. This result is an explicit verification of the asymptotic behaviour derived in Eq. (16). Eq. (28) also provides the contribution to ( ) 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 ( ) due to the reaction term, ( )/ , in the long-time limit based on Eqs (24)- (27). In Skellam's model with Dirichlet boundary condition, we obtain In the saturated growth model, by retaining the next-order term in the long-time limit, we obtain: where J ( ) = 2 It can also be checked explicitly from Eqs (24)-(27) that the reaction-dependent expressions found above by calculating ( )/ in the limit → 0, correspond exactly to − (0, ) / (0, ) , as predicted by Eq. (18). We see that in Skellam's model, the slope of the linear relationship ( ) continues to increase with time (in absolute value). In contrast, in the saturated growth model, the slope of the linear relationship ( ) decreases with time (in absolute value), and converges to as confirmed by our numerical simulations (Figures 12, 14). Figure 14 shows that the front speed law ( ) in the two-dimensional Fisher-KPP model also converges to Eq. (29) in the long-time limit, which is consistent with the fact that the saturated growth model is the linearisation of the Fisher-KPP model about the steady state = 1. At early times, however, the Fisher-KPP model is similar to Skellam's model. The slope of ( ) therefore first increases in time, before converging to that of Eq. (29) in all pore shapes (Figures 3,5-7).
Strong Allee effect. Finally, we consider an example of bistable reaction term where both the rest state = 0 and excited state = 1 are stable steady states. A classic example is provided by population growth with strong Allee effect, modelled by ( ) = ( − )(1 − ) with 0 < < 1. This model supports travelling waves with speed 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 on the shape of the initial condition [4,41]. Numerical simulations of this model are shown in Figure 13 with = 2/5 and = 1/(1 − ) = 5/3 so that ( ) ∼ 1 − as → 1 and > 0. Here too, the relationship ( ) converges to that of the saturated growth model in Eq. (29) as → ∞.
Our numerical observation that the linear relationship ( ) 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 Eq. (18), the slope of the relationship ( ) depends on (0, ). This quantity could in principle depend on how the solution profile approaches steady-state from the initial condition, where the reaction term has different contributions in the three models.

Conclusions 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 with a speed approaching the asymptotic limit Sat ( ) in Eq. (29) (blue). In these simulations, diffusivity is increased by a factor 10 compared to previous inward moving front simulations, to prevent the reaction term from dominating the evolution in the dip-filling phase.  1) is also shown for reference (green). Diffusivity is the same as in the saturated growth model for comparison in Figure 14. 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 − to the eikonal equation derived from singular perturbation theory in the low-diffusion limit [1,4,6] is in fact an exact, nonperturbative contribution originating from diffusion occurring transversally, i.e., along the wave front. For fronts moving inward into an empty domain, we identify three phases of the solution: an establishment phase, strongly influenced by conditions on the domain boundary, a travelling-wave-like phase, well described by the eikonal equation = − predicted by singular perturbation theory in the low-diffusion limit, and 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 − originating from diffusion occurring normally to wave fronts, and a time-dependent contribution originating from nonconservative changes to the solution described by the reaction term. The latter contribution becomes timeindependent if = 0 for any choice of the reaction term ( ). These behaviours are confirmed by the explicit long-time solutions of Skellam's model ( ( ) = ), the saturated growth model ( ( ) = (1 − )), as well as the strong Allee source term ( ( ) = ( − )(1 − )), 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 ( ), 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 upon that arises from transverse diffusion, and the fact that mean curvature flows evolve interfaces into shrinking circles [18][19][20][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 reactiondiffusion waves may also evolve on the surface of curved domains [50,52]. In these cases, the curvature of the surface also influences travelling wave speeds [2,51]. It is unclear how our results apply to these situations. Generalising the methods presented in this paper 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 ( ) = (1 − 2 ) [46], where the dip-filling phase is absent. Finally, generalisations of our analyses to heterogeneous excitable media [47], and to Fisher-Stefan models representing reaction-diffusion systems coupled with moving boundary conditions [48,49] are also of interest, both in situations supporting travelling-wave solutions as well as in dip-filling situations.