The goal of this blog is to understand how Melnikov theory can be used to determine the wavespeed of travelling waves in a perturbed reaction-diffusion system. In particular the focus of this blog will be on the paper Wavespeed in reaction-diffusion systems, with application to chemotaxis and population pressure (Balasuriya and Gottwald 2010). We wish to investigate perturbing a standard form of a reaction-diffusion partial differential equation (PDE) which allows us to use tools from Melnikov method to find the wavespeed of the solution. The application is again motivated by biology to describe chemotaxis of cells.
Melnikov method
Melnikov method is a technique to find chaos in dynamical systems exhibiting periodic perturbation. By constructing a “Melnikov function” we are able to predict whether a system is regular or chaotic. To determine the wavespeed of our reaction-diffusion PDE we add a small perturbation to it which allow us apply the Melnikov theory to gain an analytical solution for the wavespeed.
Reaction diffusion with strong Allee effect
Firstly, let’s explore the PDE investigated in the paper. Reaction-diffusion equations are a class of PDE that are widely used to model many physical phenomena in biology, physics, chemistry and other fields. For example, one application is in biology to model the population density of bacterium \(u(x,t)\) evolving over space and time. One particular example of a non-dimensional reaction-diffusion system in one dimension is by
\[\begin{equation} u_t = D u_{xx} + G(u) + \epsilon h(u, u_t, u_{xx}) \tag{1} \end{equation}\]
where \(\epsilon\) is a small positive number, and \(D\) is a constant diffusion coefficient. This PDE will be used to model chemotaxis of cells in this case. In Example 1 of (Balasuriya and Gottwald 2010) \(G(u) = -u(1-u)(u-\alpha)\) and \(h(u,u_x,u_{xx}) = -u^n u_x\) was defined. Hence,
\[ \begin{equation} u_t = D u_{xx} -u(1-u)(u-\alpha) -\epsilon u^n u_x, \tag{2} \end{equation} \]
where \(n\) is a positive integer and \(\alpha \in (0,1)\). We now need to define the initial and boundary conditions to find a solution for our PDE. In this case, the boundary condition is \(u(\eta) \to 0\) as \(\eta \to -\infty\) and \(u(\eta) \to 1\) as \(\eta \to \infty\). Hence, the Heaviside function \(\mathcal{H}(x)\) is an appropriate choice for the initial boundary condition.
In this example, the \(G(u)\) is a strong Allee effect because the \(\alpha\) parameter determines the split between the equilibrium occurring at \(u = 0\) or \(u = 1\). Note by equilibrium it is meant when the spatial derivatives are zero. A sketch of the strong Allee effect can be seen below.
Strong Allee effect displaying equilibrium point at 0 and \(k=1\) split by \(\alpha\).
We also note that the choice of the function \(G(u)\) excludes the simple Fisher KPP equation because Equation (1) would exhibits a continuous range of wavespeed which makes no sense to determine the wavespeed of the unperturbed scenario. In other words, we wish to solve systems of PDEs of this form that has unique wavespeeds as given in Equation (2).
Numerical solution
We first attempted to solve Equation (2) numerically using the method of lines. Method of lines is a finite difference method that discretises all variables while leaving time continuous. This results in a system of ordinary differential equation which we can solve using MATLAB’s ode solvers for example. An example of the discretisation is given below.
\[ \underbrace{\begin{bmatrix} 1 \quad 0 \quad \cdots \quad 0 \quad 0 \\ 0 \quad 1 \quad \cdots \quad 0 \quad 0 \\ \vdots \quad \vdots \quad \ddots \quad \vdots \quad \vdots \\ 0 \quad 0 \quad \cdots \quad 1 \quad 0 \\ 0 \quad 0 \quad \cdots \quad 0 \quad 1 \end{bmatrix}}_{\text{Left Hand Side (LHS)}} \begin{bmatrix} \frac{d u_1}{dt} \\ \frac{d u_2}{dt} \\ \vdots \\ \frac{d u_{M-1}}{dt} \\ \frac{d u_{M}}{dt} \end{bmatrix} = \underbrace{\begin{bmatrix} 0 \\ \vdots \\ Discretisation \\ \vdots \\ 0 \end{bmatrix}}_{\text{Righ Hand Side (RHS)}} \]
Below is the numerical solution exhibiting a travelling wave.
Solution of \(u(x,t)\) from Equation (2) with \(n = 1\), \(\alpha=0.3\) and \(\epsilon = 0.1\). The initial condition used was the Heaviside function \(\mathcal{H}(x)\) satisfying the boundary conditions \(u(x,t) \to 0\) as \(x \to -\infty\) and \(u(x,t)\to 1\) as \(x \to \infty\). We observe a distinct travelling wave from \(t=4\) onwards at a constant speed.
Analytical solution
We can now find the wave length from the numerical solution. However, it is always best to find analytical solutions to PDE whenever possible. Let’s now investigate finding wave speed using Melnikov method. The first step is to rewrite the PDE in term of wave speed. That is \(u(x,t) = u(\eta)\) where \(\eta = x-ct\) in which \(c\) is the wavespeed. Now let’s define \(w(\eta) = u'(\eta)\) and we can rewrite Equation (1) as
\[ \begin{equation} -c w = Dw' + G(u) + \epsilon h(u,w,w'). \tag{3} \end{equation} \]
The unperturbed case of Equation (1) is linked to Equation (3) by a heteroclinic connection between \((u,w) = (0,0)\) and \((1,0)\). Hence, if a nearby heteroclinic orbit persists for small \(|\epsilon|\), the perturbed stable and unstable manifolds must continue to coincide. This will allow us to using the Melnikov tools from dynamical systems to find the wavespeed of the PDE. The distance between the perturbed stable and unstable manifold need to be zero for the a heteroclinic orbit to persist. This distance is connected to the Melnikov function. To begin the analysis to find the close form of the wavespeed we need to rewrite Equation (3) as a linear system. Then utilising the following Melnikov function,
\[ M(\eta) = \int_{-\infty}^\infty \exp\left[\int_r^\eta (\boldsymbol{\nabla} \cdot \boldsymbol{f}(\boldsymbol{\bar{x}}(s))\, ds \right] (\boldsymbol{f} \wedge \boldsymbol{g}) (\boldsymbol{\bar{x}}(r))\, dr, \]
we can find an expression for the wavespeed that is given in detail in Balasuriya/Gottwald paper. Hence, the resultant wavespeed expression is given by
\[ c = \sqrt{2D} \left(\alpha - \frac{1}{2}\right) + \epsilon \left(\frac{1+2\alpha}{4}\right) + \mathcal{O}(\epsilon^2). \]
The plot of the wavespeed from the analytical and numerical solution of Equation (2) is given in the figure below.
Verification of wavespeed by varying \(D\), with \(n=1\), \(\alpha = 0.3\), and \(\epsilon = 0.1\). The dashed line is the wavespeed \(c_0\) of the unperturbed Equation (2), the solid curve is the analytical wavespeed, and the circles denote numerically computed wavespeed calculated from a simulation of the full partial differential equation.
From the figure above we can see that the numerical and analytical solutions agree. Hence, the Melnikov method is a reliable technique in finding the wavespeed of reaction-diffusion equations under the condition that wavespeeds are unique.