Front stability of infinitely steep travelling waves in population biology

TL;DR

Full 2D time–dependent level–set numerical solutions indicate that invasive travelling waves are stable to small amplitude lateral perturbations, whereas receding travelling waves are unstable.

Abstract

Reaction-diffusion models are often used to describe biological invasion, where populations of individuals that undergo random motility and proliferation lead to moving fronts. Many models of biological invasion are extensions of the Fisher-KPP model that describes the evolution of a 1D population density as a result of linear diffusion and logistic growth. In 2020 Fadai introduced a new model of biological invasion that was formulated as a moving boundary problem with a nonlinear degenerate diffusive flux. Fadai's model leads to travelling wave solutions with infinitely steep, well-defined fronts at the moving boundary, and the model has the mathematical advantage of being analytically tractable in certain parameter limits. We aim to provide general insight by first presenting two key extensions by considering: (i) generalised nonlinear degenerate diffusion with flux; and, (ii) solutions describing both biological invasion, and biological recession. We establish the existence of travelling wave solutions for these two extensions, and then consider stability of the travelling wave solutions by introducing a lateral perturbation of the travelling wavefront. Full 2D time-dependent level-set numerical solutions indicate that invasive travelling waves are stable to small lateral perturbations, whereas receding travelling waves are unstable. These preliminary numerical observations are corroborated through a linear stability analysis that gives more formal insight into short time growth/decay of wavefront perturbation amplitude.

Figures (8)

• Figure 1: Schematic travelling wave solutions: (a) Fisher--KPP model on $-\infty < x < \infty$; (b) Porous--Fisher model on $-\infty < x < \infty$; (c) Fisher--Stefan model on $x < L(t)$ ; and (d) the Porous--Fisher--Stefan (PFS) model on $x < L(t)$. Arrows on the travelling wave profiles in (a)--(d) indicate that the Fisher--KPP and Porous--Fisher models describe biological invasion only, whereas the Fisher--Stefan and PFS models can describe biological invasion and biological recession. Schematics in (e)--(f) illustrate four possibilities for the linear stability of travelling wave solutions of the PFS model with $u \to 1$ as $x \to -\infty$. The profile at $t=0$ (red) shows the location of the moving boundary where $u=0$ including a small amplitude transverse perturbation. Profiles at $t_1 > 0$ (green) and $t_2 > t_1 > 0$ (blue) and shows the evolution of the moving boundary. Arrows in (e)--(h) show the direction of increasing time. Profiles in (e)--(f) are stable invading and receding fronts, respectively. Profiles in (g)--(h) are unstable invading and receding fronts, respectively.
• Figure 2: Numerical solution of Equations (\ref{['Eq:PFS']})--(\ref{['Eq:PFS MovingBC']} on a rectangular domain with $-10 < x < 10$ and $0 < y < 10$, and $u(x,y,0) = 1$ for $x < 0$ with $\partial \Omega(0)$ along the vertical line at $x=0$. Numerical simulations are performed with periodic boundary conditions along the horizontal boundaries at $y=0$ and $y=10$, and we set $u=1$ along the vertical boundary at $x=-10$. Numerical solutions are obtained on a $401 \times 201$ uniform finite difference mesh. In each plot we show solutions at $t=0, 30, 60$ and $90$ with the arrow showing the direction of increasing time. Results in (a)--(c) correspond to $\kappa=0.1$ whereas results in (d)--(f) correspond to $\kappa = -0.1$. The value of the exponent $m$ is indicated on each column. We report the value of the front velocity at the end of the simulation, $c$.
• Figure 3: Phase planes and associated trajectories for Equations (\ref{['eq:theta_prime']})--(\ref{['eq:psi_prime']}). Phase planes in (a)--(c) correspond to $m=0.5, 1.0, 2.0$ as indicated, and each phase plane shows the location of the equilibrium points at $(0,0)$ and $(1,0)$ with black discs. Each phase plane shows trajectories associated with travelling wave solutions for $c=-0.5, -0.25, 0, 0.25, 0.5$ in purple, yellow, red, green and blue, respectively. Each trajectory intersects the $\psi$ axis at a particular point shown in a coloured dot. The coordinates of this point allow us to calculate $\kappa$ according to $\kappa = -c(1+m)/ \psi(0)$. For $c=-0.5, -0.25, 0, 0.25, 0.5$ and $m=0.5$ we obtain $\kappa=-0.59,-0.38,0.0,0.79, 2.82$, respectively. For $c=-0.5, -0.25, 0, 0.25, 0.5$ and $m=1.0$ we obtain $\kappa=-0.63,-0.42,0.0,1.05, 5.4$, respectively. For $c=-0.5, -0.25, 0, 0.25, 0.5$ and $m=2.0$ we obtain $\kappa=-0.69,-0.49,0.0,1.87, 290241.59$, respectively.
• Figure 4: Relationships between $c$, $\kappa$ and $m$ deriving using the phase plane. Results in (a) show $c$ as a function of $\kappa$ for $m=0.5, 1.0, 2.0$ in green, blue and red curves, respectively, and the arrow shows the direction of increasing $m$. Results in (a) are obtained by solving Equations (\ref{['eq:theta_prime']})--(\ref{['eq:psi_prime']}) for certain values of $c$ and $m$, and then applying the boundary condition $\kappa =-c(1+m)/ \psi^*$ to estimate $\kappa$. Results in (b) are for the special case $m=1$ where the solid blue curve is identical to the solid blue curve in (a), and the dashed red curve is the perturbation result, Equation (\ref{['eq:c_kappa_perturbed_relation']}). The inset shows a plot that focuses on the dashed rectangular region near $c=0$ where we see that the perturbation result is visually indistinguishable from the numerically--generated phase plane results at this scale. Each curve in (a)--(b) obtained using the phase plane is constructed by using 50 equally--spaced values of $c$ for each $m$, and then using the phase plane to solve for $\kappa$. Interpolating the $(c, \kappa)$ values gives the continuous curves in (a)--(b).
• Figure 5: Numerical exploration of front stability with $\kappa = 0.1$. Profiles in (a)--(c) show three perturbed travelling wave profiles for $m=0.5, 1.0, 2.0$, respectively. Each perturbed travelling wave profile is computed by first solving the $\mathcal{O}(1)$ boundary value problem to give the shape of the unperturbed travelling wave that is then shifted so that $\partial \Omega(0)$ coincides with the vertical line at $x=5$ before a perturbation with wavenumber $q=4\pi/5$ and amplitude $\varepsilon = 0.1$ is added to give the sinusoidal profiles in (a)--(c) at $t=0$. Results in (d)--(f) show the solution of Equations (\ref{['Eq:PFSphi']})--(\ref{['Eq:PFSphi MovingBC']}) at $t=50$ where, in each case, we can see that each front moves in the positive $x$--direction and amplitude of the transverse perturbation clearly decays. Plots in (g)--(i) show the evolution of $\partial \Omega(t)$ at $t=0, 2, 4, 6, 8$ with the arrow showing the increasing value of $t$. Note that the plots in (g)--(i) focus on the interval $4 \le x \le 6$ so that the details of the interface evolution are clear. Comparing these plots of $\partial \Omega(t)$ suggests that these travelling waves with $\kappa = 0.1$ are stable as the perturbation amplitude appears to decay with time. All numerical calculations are performed on a $10 \times 10$ domain uniformly discretised with a $201 \times 201$ mesh.