May 20, 2018, Sunday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Propagation of fronts in the Fisher-Kolmogorov equation

From MathBio

Jump to: navigation, search


This page details a model governing behavior of fronts in heterogeneous environments, as given in paper Propagation of fronts in the Fisher-Kolmogorov equation with spatially varying diffusion by Curtis and Bortz[1]. Paper provides an analytic framework that describes propagation of the front up to the minimum of the diffusion coefficient and shows that standard traveling coordinate frames do not properly describe front propagation.


Wave travelling in positive direction
Reaction-diffusion systems are systems involving components locally transformed into each other by chemical reactions and transported in space by diffusion. They arise, quite naturally, in chemistry and chemical engineering but also serve as a reference for the study of a wide range of phenomena encountered beyond the strict realm of chemical science such as environmental and life sciences [2]. Mathematically, reaction–diffusion systems take the form of semi-linear parabolic partial differential equations. They can be represented in the general form

where each component of the vector u(x,t) represents the concentration of one substance, D is a diagonal matrix of diffusion coefficients, and f(u) accounts for all local reactions. The solutions of reaction–diffusion equations display a wide range of behaviours, including the formation of travelling waves and wave-like phenomena as well as other self-organized patterns like stripes, hexagons or more intricate structure like dissipative solitions. [3] The choice f(u)=u(1-u) yields Fisher's equation . Fisher proposed this equation to describe the spatial spread of an advantageous allele and explored its travelling wave, or propagating front, solutions. [4]


Front propagation in reaction-diffusion equations (RDEs) is an important topic in the physical and biological sciences. Interest in the behavior of fronts in heterogeneous environments has increased over the last several years. In particular, traveling waves, or propagating fronts, through certain classes of time and spatially varying environments have been well studied. Nonlinear, density-dependent diffusion, where the diffusion coefficient depends on u, has also been thoroughly examined. In the biological literature, spatially heterogeneous diffusion coefficients a(x) are discussed repeatedly, but continuously varying models are rarely investigated rigorously. A morphogenesis phenomenon in which a(x) has been given serious consideration is in a general RDE model for the slug stage of Dictyostelium discoideum. The resulting diffusion profile across a gap junction is a hyperbolic cosine, i.e. a(x)=cosh(x). Moreover, recent work has considered the impact of letting a(x)=D+\eta x^{2} (for D constant and η small) in the context of the avascular growth phase of cancer as well as on Turing bifurcations of standing wave solutions to RDEs.

Mathematical model

This paper considers the case of a heterogeneous Fisher-Kolmogorov (FK), or Kolmogorov-Petrovskii-Piskunov (KPP), equation where heterogeneities are represented by spatially varying diffusion, i.e.


And f(u)=u(1-u).

The form of a(x) discussed in this page is a(x)=x^{2}+\epsilon . This choice offers a good canonical model for a diffusion-mediated barrier. For example, one could use this mathematical structure to model a geographic barrier such as a mountain. Further, this choice also allows us to study the effect of strongly varying diffusion.

Initial condition taken is the step initial condition

u(x,0)=\left\{{\begin{array}{cc}1,&x\leq x_{c}(0)\\0,&x>x_{c}(0)\end{array}}\right.

Where x_{c}(t) denotes the location of the front at time t.

Analysis from the paper

The Soft Front Approximation

When, we assume that the solution has a soft front, i.e., a condition defined via the asymptotic relationship

a(x)u_{{xx}}\ll a'(x)u_{x}

This asymptotic condition is denoted as the Soft Front Approximation (SFA). Using the SFA, we now must solve the semilinear hyperbolic equation,

u_{t}=a(x)u_{{xx}}(x)+a'(x)u_{x}+f(u)\thicksim a'(x)u_{x}+u(1-u)

Then using method of characteristics, we get the solution

u(x,t)\sim {\frac  {u(xe^{{2(t-t_{0})}},t_{0})}{1+u(xe^{{2(t-t_{0})}},t_{0})(e^{{t-t_{0}}}-1)}}\qquad (1)

which follows the characteristic curve x(t)=x_{0}e^{{-2(t-t_{0})}}. In the case that x_{0}<0, we see that information propagates towards the origin as desired for a traveling front. However, if x_{0}>0, then all information again propagates to the origin, and thus the soft front model cannot describe a front propagating past the origin. Hence, we have a mechanism that explains front propagation that follows a decreasing diffusion coefficient; i.e., a(x) is strictly decreasing when x_{c}(t)<0. As can be seen from (1), the front does not propagate as a traveling wave; i.e., there is not a constant wave speed such that u(x,t)\sim u(x-ct). However, it is not hard to see that the front travels along curves of the form

\ln {(x(t))}+2(t-t_{0})=\ln {(x_{c}(t_{0}))}\qquad (2)

In this log-transformed spatial coordinate we can see the front propagates with speed c=2. We refer to (2) as a traveling wave coordinate (TWC), and we see that the TWC and the SFA allow us to compute a generalized notion of wave speed. They argue that the SFA holds in a region around x_{c}(t) and this region must be matched, via intermediate layers, to the far field of u(x,t), which should have a much steeper decay profile. To study this regime, they linearize the FK equation around u=0; i.e., let u={\tilde  \epsilon }v and then collect all terms in {\tilde  \epsilon }, so we get


Assuming that v is given by the WKB ansatz v(x,t)=A(x,t)e^{{\phi (x,t)}} we get the leading order problem

\phi _{t}+a(x)\phi _{x}^{2}+1=0

Using the method of characteristics and noting that above equation is a Hamilton-Jacobi equation, with Hamiltonian H(p,x)=1+a(x)p^{2}, we get that characteristic associated with the choice of diffusion coefficient a(x)=x^{2}+\epsilon are given by

x_{0}={\frac  {1}{2}}(x+{\sqrt  {x^{2}+\epsilon }})e^{{\mp 2{\tilde  H}t}}-{\frac  {\epsilon }{x+{\sqrt  {x^{2}+\epsilon }}}}e^{{\pm 2{\tilde  H}t}}

Thus we see that for |x|\gg {\sqrt  {\epsilon }}, the characteristics are to leading order given byFailed to parse (lexing error): x_0∼xe^{\mp 2\tilde H t } .

In the FK equation, the minimum of the diffusion coefficient represents a turning point; i.e., the sign of a a'(x) changes through the minimum. Through the turning point, or when |x|\ll {\sqrt  {\epsilon }}, we get the leading order behavior

x_{0}\sim \mp {\sqrt  {\epsilon }}\sinh {(2{\tilde  H}t)}+(x+{\frac  {x^{2}}{2{\sqrt  {\epsilon }}}})\cosh {(2{\tilde  H}t)}

The WKB analysis shows how the SFA eventually breaks down and as the front approaches the origin, something like a shock, or steepened front, must form.

Softening Sharp Initial Contditions

We assume throughout this section that x_{c}(0)\gg 1. The asymptotic condition describing the sharp front is that

a(x)u_{{xx}}(x)\gg a'(x)u_{x}

We also choose a parameterL\gg x_{c}\gg 1 so that u_{x}(\pm L,t)=0. We choose Neumann boundary conditions to allow for analytical tractability. Introduce the fast time T={\frac  {t}{\epsilon }} and using the ansatz,

u=u_{0}(x,T,t)+\epsilon u_{1}(x,T,t)\cdots

we get the equations

\partial _{T}u_{0}=\partial _{x}(a(x)\partial _{x}u_{0})
\partial _{T}u_{1}=\partial _{x}(a(x)\partial _{x}u_{1})+f(u_{0})-\partial _{t}u_{0}

Leading order equation can be solved using separation of variables, which leads to expansion for u_{0}

u_{0}(x,t,T)=\sum _{{n=0}}^{{\infty }}{\sigma _{n}(t)\phi _{n}(x)e^{{\lambda _{n}T}}}

where \phi _{n} and \lambda _{n} solve the Sturm-Liouville problem

\partial _{x}(a(x)\partial _{x}\phi _{n}(x))=\lambda _{n}\phi _{n},\quad \partial _{x}\phi _{n}(\pm L)

with initial conditions \lambda _{0}=0,\quad \phi _{0}={\sqrt  {1/2L}},\quad \sigma _{0}(0)={\frac  {x_{c}+L}{{\sqrt  {2L}}}}\quad \sigma _{n}(0)={\frac  {a(x_{c})}{\lambda _{n}}}\partial _{x}\phi _{n}(x_{c}/\epsilon ) Sturm-Liouville theory ensures the functions \phi _{n} form a complete, orthonormal set with respect to L_{2} norm. Since we assume that a(x)>0, it easily follows that \lambda _{n}\leq 0

Moving to the second term u_{1}(x,T,t), we see, using Duhamel’s principle, that we can write u_{1} as

u_{1}(x,T,t)=\int _{{0}}^{{T}}{\sum _{{n=0}}^{{\infty }}{\gamma _{n}(t,s)\phi _{n}(x)e^{{-|\lambda _{n}|(T-s)}}}}


\gamma _{n}(t,s)=\int _{{-L}}^{{L}}{(-\partial _{t}u_{0}+u_{0}(1-u_{0}))\phi _{n}dx}

Since \lambda _{0}=0, one can see that a possibility for a secularity, which means the asymptotic series becomes invalid on O(1) time scales, arises from computing

\int _{{0}}^{{T}}{\lambda _{0}(t,s)ds}

Thus enforcing below condition removes secularity

{\dot  \sigma }_{0}=\sigma _{0}(1-\phi _{0}\sigma _{0})

which has the solution

\sigma _{0}(t)={\frac  {1}{\phi _{0}+(1/\sigma _{0}-\phi _{0})e^{{-t}}}}

\sigma _{n}(t) can also be calculated:

\sigma _{n}(t)={\frac  {\sigma _{n}(0)e^{t}}{(1+\sigma _{0}(0)\phi _{0}(e^{t}-1))^{2}}}

Thus, on timescales T=O(1/\epsilon ), or t=O(1), one has

u(x,T,t)=u_{0}(x,T,t)+O(\epsilon )

From this, since for n\geq 1 one has by orthonormality

\int _{{-L}}^{{L}}{\phi _{n}(x)dx}=0

it follows that the average of u(x,T,t) to leading order is given by \sigma _{0}(t)\phi _{0}, or

{\frac  {1}{2L}}\int _{{-L}}^{{L}}{u(\xi ,T,t)d\xi }=<u>\sim \sigma _{0}(t)\phi _{0}={\frac  {1}{1+{\frac  {L-x_{c}}{L+x_{c}}}e^{{-t}}}}

This result shows that the average should increase exponentially fast on time scales of O(1) or less, which amounts to showing that the transition region of the solution u becomes smoother and “softer” rapidly. Having done this, the solution then enters the SFA regime.

Front Beyond The Turning Point: Reduction To A Stationary Equation

Since the SFA only leads to a model that describes propagation of a front up to the turning point, we must now find some other means of trying to describe propagation of the front past the turning point. To do this, by taking x\gg 1, so that a(x)\thicksim x^{2}, we suppose that u=u(\eta (x,t)), which means the FK equation becomes

\eta _{t}{\frac  {du}{d\eta }}=(x\eta _{x})^{2}{\frac  {d^{2}u}{d\eta ^{2}}}+(x^{2}\eta _{x})_{x}{\frac  {du}{d\eta }}+f(u)

By choosing

\eta _{x}=\pm {\frac  {1}{x}},\quad \eta _{t}=c

we get the stationary FK equation

{\frac  {d^{2}u}{d\eta ^{2}}}+(-c\pm 1){\frac  {du}{d\eta }}+f(u)=0

with coordinate

\eta (x,t)=\pm \ln {|x|}+ct

The variable \eta is another instance of a TWC, and so again we see that strongly varying diffusion requires a generalization of the definition of a traveling front.

Numerical results

Figures 1 and 2 depict the development of the solution in the (t, x) plane as a surface and as a contour plot, respectively. Figure 3 shows the front before, at, and after the turning point. In Figure 4, we compare the numerical propagation of the front and our asymptotic the two curves represent the function x_{{c,n(t)}} at which u_{{num}}(x_{{c,n(t)}},t)=1/2, where u_{{num}}(x,t) denotes the numerical approximation to u(x,t). In the case of the asymptotic curve, we take t_{0}=0, and use the same step initial condition as used in the numerics in (4). The agreement between the SFA and numerics is convincing, and thus we have confirmation that the SFA is valid up to the turning point. As can be seen the SFA breaks down at the turning point. From Figures 1, 2, and 3, we see that the front steepens substantially.

Picture1.png Picture2.png

Pic3.png Pic4.png

In Figure 5, we plot the time duration that the simulated front spends near x=0 as a function of 1/\epsilon , where we define near as |x|<0.4. Note the value 0.4 is the width of the spatial mesh step in the numerics, and thus the smallest scale on which phenomena can be distinguished. Figure 5 shows that a least squares fit of the times to a curve growing like {\frac  {1}{{\sqrt  {\epsilon }}}} is accurate.


On the other side of the turning point, as seen in Figure 2, the solution propagates in a manner consistent with the log-transformed TWC, i.e., information is propagating along curves of the form \eta =-\ln {|x|}+ct


After diffusion smooths the initial conditions such that the SFA is valid, the choice of a strongly varying diffusion coefficient with a global minimum implies the following

  • The definition of a traveling front must be generalized via the Traveling Wave Coordinate
  • The minimum, or turning point, of a(x) causes the formation of shocklike behavior and leads to a trapping of the front on asymptotically long time scales.
  • The behavior of the front on either side of the turning point is fundamentally different, and on either side, the TWC is necessary to describe dynamics.
  • After the front is past the turning point, the TWC allows us to transform into a stationary FK equation, thus simplifying all subsequent analysis


Mountain barrier
  • Quadratic spatial diffusion could be used to model barriers such as mountains, and results will give estimates of a time to barrier transit. For example, the invasion of the Midwestern United States by gypsy moths has been extensively studied[5]. It is also well known that the moths’ spatial diffusion rate is tightly correlated with the local habitat. By estimating the population growth rates and fitting a quadratic curve to the local diffusion rates from capture-mark-recapture experiments, the method developed in this paper could be used to estimate the waiting time until the moths reach a certain habitat.
  • Nontrivial spatially varying FK equation, studied here, could be reduced to a much simpler stationary problem. This will allow for a greater degree of control and flexibility in curve-fitting routines as it will only require a snapshot at a single point in time. With the gypsy moth invasion example, one could envision estimating the wave speed from a single snapshot.
  • It would be of particular interest to develop asymptotic matching techniques that would allow for connecting the various regimes of the front, i.e., the sharp to soft front transition, and then the trapping layer around the minimum of a(x). Likewise, it would also be of interest to develop higher order expansions for the speedcin the SFA regime.


  1. C. W. Curtis and D. M. Bortz, Propagation of fronts in the Fisher-Kolmogorov equation with quadratically varying diffusion. Phys. Rev. E, 86(6), 066108, Dec. 2012.
  4. Fisher, R. A., The genetical theory of natural selection. Oxford University Press, 1930. Oxford University Press, USA, New Ed edition, 2000
  5. P. C. Tobin, S. L. Whitmire, D. M. Johnson, O. N. Bjørnstad, and A. M. Liebhold,Ecol. Lett.10, 36 (2007)