More than Skew: Asymmetric Wave Propagation in a Reaction-Diffusion-Convection System

Convection-induced instability in reaction-diffusion systems produces complicated patterns of oscillations behind propagating wavefronts. We transform the system twice: into lambda-omega form, then into polar variables. We find analytical estimates for the wavefront speed which we confirm numerically. Our previous work examined a simpler system [E. H. Flach, S. Schnell, and J. Norbury, Phys. Rev. E 76, 036216 (2007)]; the onset of instability is qualitatively different in numerical solutions of this system. We modify our estimates and connect the two different behaviours. Our estimate explains how the Turing instability fits with pattern found in reaction-diffusion-convection systems. Our results can have important applications to the pattern formation analysis of biological systems.


I. INTRODUCTION
Reaction-diffusion systems have been studied since Turing suggested them as an explanation of morphogenic patterns [1].Mathematical biologists have shown that Turing-type models can mimic patterns during developmental processes, such as pigmentation patterns on animals [2], feather germ patterns form in a hexagonal array behind a propagating maturation front [3] or traveling stripes on the skin of a mutant mouse [4].More recently, this system has been found in mouse hair follicle dis-tribution [5].Here we investigate the effect of adding convection on a general model of pattern formation.Pattern formation occurs as organisms grow, and this increase in size could be modelled as convection.More generally, adding convection to the reaction-diffusion system is a disturbance to the system.This means our approach can be considered a stability analysis of the original system.We see that the convection gives rise to Turing-type pattern more readily.That the system produces pattern more easily under a disturbance gives us more confidence in the mechanism for biological applications.
The reactions often chosen for study are those proposed by Schnakenberg [6], since he demonstrated that they originated from the simplest two-species chemical reaction capable of producing limit cycles in an ODE.Cubic auto catalysis is a model for several chemical reactions [7].As discussed by Schnakenberg, the limit cycle requires an unstable spiral and bounding on the instability.These two properties give functions suitable for pattern formation in a PDE.
More recently, Rovinsky and Metzinger [8] showed that pattern could be induced by convective flow.They employ the more complicated, but chemically valid, chlorite-iodide-malonic acid system.Subsequently, diffusion was introduced to this system and instability observed experimentally [9].Satnoianu et al. [10] showed Figure 1.Pattern found for a diffusion system with convection and limit-cycle reaction kinetics (1).The initial disturbance propagates and becomes pronounced, forming a regular pattern with aligned oscillations.The propagation is linear, forming a V-shape.This is a numerical solution using NAG D03PCF, plotting species u with µ = 1.1, D1 = 0.01, D2 = 0.01, ρ = 1.The reactants are initially at steady state: (u, v) = (1/µ, µ), with a small disturbance at x = 0.The boundaries are held at zero derivative: ux, vx = 0.
that an initial point disturbance in this system can cause a pattern to propagate.The pattern is an oscillatory form, two wave packets joined together, and increasing in width.We show this pattern in Figure 1.This is the propagation of a small disturbance from the steady state, and so the initial conditions are the steady state with a unit width disturbance an order of magnitude smaller than the pattern at x = 0 in the first species u only.This disturbance is apparent in Figure 1.For simplicity we imagine an infinite domain; for numerical solutions we keep the boundaries at the steady state value.These conditions are retained throughout this study.This propagating pattern shown in Figure 1 is interesting but poorly understood, in the general case.
The speed of propagation has been given for a twospecies system with no convection [11] and one with equal convection [12].A related theoretical study considers the relationship between the onset of the instability and the longer-term behaviour of the non-equilibrium state [13].Here the simplest reaction-diffusion system with oscillatory kinetics was considered, and complex behaviours were found.
In a previous paper, we examined a similar but simpler system [14] using λ-ω reactions because their symmetry is such that the chemical species convert readily into polar form.Murray's book Mathematical Biology gives a good introduction to λ-ω systems and their application [15].We analysed the pattern propagation as two travelling waves of oscillations emanating outwards from a central axis.Our analytical predictions were confirmed by numerical observation.This result in our previous paper led us to the technique we present here for the general analysis of a physical reaction-diffusion system with convection.We focus our attention to the physical reaction-diffusion system introduced by Satnoianu et al. [10].This is the reaction-diffusion-convection system with Schnakenberg reaction kinetics.
Our previous success with the simpler λ-ω system suggests that the same approach will be fruitful: transformation of the dependent variables into polar form.However, the Schnakenberg reactions lack the symmetry of the λ-ω functions and so a further technique is required.We employ a simple technique suggested by Murray [15, in a novel situation.This is a matrix transformation, converting the linearised Schnakenberg reactions into λ-ω form.Having found a symmetric representation, we convert into polar form, then apply a Fisher travelling wave analysis.This gives predictions for the rate of expansion of the envelope of pattern.We refer to this as the wave speed of the boundaries.
For this previous system we found two distinct wave speed classes with the switch a function of the strength of convection.In the absence of any other information, we predict that the Schnakenberg system will behave in the same way as the λ-ω system.We investigate numerically using NAG D03PCF.We start with the system at the steady state, and make a small perturbation in u at x = 0. To avoid boundary effects, we hold the derivative at zero on the boundary.We look at three values for µ, the source supply term.For each µ, we vary both ε, the difference in diffusion, and γ, the convection strength.In each simulation, we automatically determine the lines of onset of pattern in x-t space: the left and right wave speeds.
Our analytical approach gives the general estimate for the wave speeds.However, the Schnakenberg system has a different switch between behaviours.In particular, the left wave speed only has one class of behaviour: the switch found in the λ-ω system is not present.We investigate further and find an underlying behaviour consistent with the switch found previously.However, left-hand stability is retained in the Schnakenberg system.

II. SCHNAKENBERG REACTION-DIFFUSION-CONVECTION SYSTEM
The most general form of the reaction-diffusionconvection system is: with the diffusion coefficients D i positive.We define the system in abstract form to demonstrate the generality of the results that follow.We remove one of the convective terms by a simple change of coordinates from z to Z = z − c 1 τ , with T = τ : We can choose the convection coefficient to be positive with an even simpler coordinate change, Y = sign(η)Z: so ρ = sign(η)η, and certainly ρ > 0. (There is the trivial case too, c 1 = c 2 , for which we ignore this step.)We expect that a difference in size in morphogens will give a difference in diffusion coefficients D i and will represent a biologically realistic pattern formation system [5].For passive convection or advection given here by ρ, we look for a flow of liquid in the system.The morphogen v is carried by this flow.The liquid flow has no effect on the morphogen u.The last step in the reduction is to remove the first diffusion coefficient with y = Y / √ D 1 : The above applies in general to this reaction-diffusion-convection system.Now we introduce the Schnakenberg reaction scheme: Using the law of mass action, this reaction scheme translates into our reactions: with µ the constant, positive source term.The rate constant k (a third-order reaction coefficient) and the sink rate constant m (first-order rate coefficient) are always positive.In full, the system is We can remove some of these constants by rescaling our dependent and independent variables together.We choose ũ = k/m.u,ṽ = k/m.v,t = mT , x = √ my ũt = ũxx + μ − ũṽ 2 ṽt = εṽ xx − γṽ x + ũṽ 2 − ṽ where γ = γ/ √ m and μ = k/mµ/m.The rescaled reactions are simpler: with only a single source parameter, µ.From this point we drop all the tildes for clarity: with The remaining constants are: ε denoting the ratio of the diffusion coefficients, γ the difference in convection, and µ the relative strength of the chemical source.

III. TRANSFORMING THE REACTIONS INTO λ-ω FORM
We have so far simplified the differential operators in our reaction-diffusion-convection model.However, the reaction terms remain complicated.Our approach is suggested by Murray [15, when considering a related problem, that of travelling wave trains.We consider the behaviour close to the steady state and so linearise about it.This straightforward transformation yields a "normalised" form of the reaction system.
To simplify the exposition, we first consider only the ODE formed using the rescaled Schnakenberg reactions: This system is known to form limit cycles from µ < 1 down to µ 0.90032 [16].2. Phase space for the limit-cycle reaction (3), given by numerical solution.The phase curve spirals out from the steady state to meet the limit cycle (broad loop).The Schnakenberg (a) system is highly elliptical, even at the start of the trajectory, by the steady state.The transformed Schnakenberg (b) shows a more circular behaviour initially, but the limit cycle shape remains more complicated.The parameter here is µ = 0.95.The solution was found using Matlab ode15s.
We find the non-trivial steady state at (1/µ, µ).We linearise the reactions about this steady state, yielding: The Hopf bifurcation point is µ 2 = 1.We look for a transformation matrix P that will convert the Schnakenberg matrix M into the symmetric λ-ω form, so that P −1 MP gives: Employing some algebraic manipulation, we find that α must be equal to (1 − µ 2 )/2.For µ near the Hopf bifurcation, α will always be small.Since the determinants of M and L must be the same, then α 2 + β 2 = µ 2 , and so β is also determined.Since α will be small, β will be similar to µ.We find a transformation matrix: and we confirm that L = P −1 MP.This matrix P is not unique; any rotation or enlargement of it is also suitable.We can use this transformation as a change of variables, namely u = P ū.This transformation only regularises the linear part of the Schnakenberg functions.We hope the transformed system will exhibit symmetry close to the steady state.However, away from this point, the non-linear terms will have a more significant contribution.Therefore we expect increasing asymmetry as the solution moves away from the steady state.
We can see the effect this transformation has on the ODE phase plane in Figure 2. The original Schnakenberg system Figure 2(a) is highly elliptical, even at the start of the trajectory, by the steady state.The transformed Schnakenberg Figure 2(b) shows a more circular behaviour initially, but the limit cycle shape remains more complicated.This is as we expected.
We can now write the PDE in matrix notation: Now we apply the change of variables to this PDE system.This gives the differential equation as: where H represents the higher-order terms.Here we have introduced a new variable δ = ε − 1, the difference between the diffusions, to collect terms.We concern ourselves only with the neighbourhood of the steady state and so neglect these higher-order terms.From now on we are dealing with the linear approximation of the system.The next step is to pre-multiply the entire differential equation ( 4) by the inverse transformation matrix P −1 , yielding: We already know that P −1 MP = L, and the other two are calculations are fairly straightforward, yielding: We can simplify the expression of the system: This clearly shows the connection between the diffusion difference δ and the convection γ.This combination term acts as a disturbance to the simple reaction-diffusion form.Since the reaction term L is now in classic λ-ω form, any difference in predicted behaviour is probably due to the modified differential form.It is possible that the non-linear terms make some contribution to the initiation of pattern too.

IV. TRANSFORMING THE CHEMICAL SPECIES INTO POLAR FORM
The purpose of the previous section was to manipulate the system into a particular form.This λ-ω form has a substantial symmetry to it, and is very suitable for the next stage in our analysis.This next step is a further transformation of the variables.We convert the newly-transformed dependent variables (ū, v) into a new coordinate system, a polar form (r, θ) as follows: The purpose of this transformation is to reduce the system to one dependent variable, by assuming that r is the significant term, and that the effect of θ is negligible.There are some useful identities which we employ: Utilising these, we find the following: Observe that the second trigonometric term (in brackets) is the differential of the first one.

A. Simple case
We consider the trivial case δ = 0, γ = 0, and assume θ x = θ xx = 0. We assumed earlier that we could dispense with the nonlinear terms, and thus we have: Then the θ equation is separable and we find θ ≈ −βt, to within a constant.The r equation is certainly suitable for a Fisher-type propagating wave solution analysis.This yields a minimal travelling wave speed estimate of ±2 √ α.In simpler systems it has been shown that this minimal speed will be achieved.However, our more complicated system does not lend itself so readily to analysis [11].It is also possible that nonlinearity may have an effect [17].To continue our analysis, we assume that our minimal wave speed estimate will be achieved in our system.

B. Full system
We complete the transformation to the polar dependent variables: we first consider the simple case where θ x = θ xx = 0, giving: The radius (r) equation determines the onset of instability.That is, an increase in distance from the steady state is exactly an instability.We make the assumption that the radius equation is the important one.It is not certain that the radius equation drives the system behaviour.However, later on, in Section VI, we will see that this appears to be the case.Therefore, for this analysis, we treat θ as if it is a constant.We carry out the standard Fisher travelling wave analysis on this equation.We say r(x, t) = R(s) at the wavefront, where s = x−ct is the travelling wave variable.We substitute this variable into the PDE, which reduces the system to an ODE: This resultant ODE is solved by an exponential, R = e ms .We substitute the exponential for R, giving a quadratic expression in m: The value for m is found by solving this quadratic.However, we are not interested in m, but the possible values for c.
The travelling wave argument is that if m is complex, then we have a spiral and so r must become negative at some point, which is unrealistic.Therefore m must be real.Since m is real, from the solution of the quadratic: In the case that the right hand side is negative, then the inequality is always satisfied: m is real regardless of the choice of c.
The next claim is that the minimal allowable wave speed c is achieved.That is, we choose c to equalise the condition.This then gives us: where c is the estimate for the wave speed.In the case where the right hand side of the inequality is negative, (1 + δa cos θ)α < 0, then the closest to the minimum is This special form is only a function of γ, although ε and µ are involved in the condition.To achieve a spreading pattern, in this case we certainly need to find different angles θ for left and right wave speeds.
The expression a cos θ can be re-written as [1 + p cos (2θ − q)] /2 where p 2 = 1 + [(1 − α)/β] 2 and tan q = (1 − α)/β.The magnitude p simplifies to √ 2µ/β.For values of µ near the Hopf bifurcation point, µ = 1, β ≈ µ and so p ≈ √ 2 ≈ 1.41.The value of θ is not important at present, so we can work with φ = 2θ − q instead.The wave speed expression then becomes For equal diffusion on both species δ = 0, the wave speed estimate is purely dependent on the convection term.Similarly, for no convection difference γ = 0, the wave speed estimate is only diffusion-dependent.Therefore, when convection is small we see propagation caused primarily by diffusion, and pattern primarily as the solution of the temporal ODE [18].This gives an oscillation at the onset of pattern, which is equivalent to θ increasing.Our small parameter estimate is for average values of the trigonometric functions: mean(cos φ) = 0.This assumes that φ and so θ is a linear function of t, which is the first solution we found.We now revert to using ε for simplicity.This reduction gives a relatively simple equation, but one which combines all the parameters in the system: This averaged wave speed gives no spread in the special c case.
For large convection difference γ 1, we consider the extremal values of the trigonometric function: cos φ = ±1.This situation yields the following wave speeds: V. NUMERICAL EXPERIMENTS We have generated an expression for a wave speed estimate, and explored this estimate briefly.We now conduct some numerical experiments to validate the estimate.There are three parameters, each of which need to be varied independently of the others.This requires an extensive set of experiments.We start by choosing a particular value of µ, in this case µ = 0.95, giving α = 1 − µ 2 /2 ≈ 0.05 (and then β ≈ 0.949, p ≈ 1.416).
This value is in the unstable region for the ODE, and so pattern should form readily.We vary ε and γ over a wide range and measure the incident wave speeds in each case.We see if the actual travelling wave speeds found corresponds to our estimated speeds (11), (12).
The results are consistent with the general analytical estimates (10), as we can see in Figure 3.However, the behaviour is not as we predicted in the specific estimates.For the previously-analysed λ-ω system we saw a transition from low convection to high convection [14].In this Schnakenberg system we see no such transition for the left wave speed: the system remains at the low-parameter estimate.In contrast, the right wavefront almost immediately diverges from our lowparameter estimate.The right wave speed approaches In contrast, the right wave speed is far from the low parameter estimate, and fairly close to the high parameter estimate.While the behaviour corresponds to the analysis, the behaviour of θ is different to expectations.In the left wavefront speed, the darker surface is the low parameter estimate: γ/2 − √ 1 + ε 1 − µ 2 ; The lighter surface is the high parameter estimate: In the right wavefront speed, the lighter surface is the high parameter estimate: The similar dark surface is a less extremal estimate, with p cos φ = 1, giving the surface as γ + √ 2ε 1 − µ 2 ; the mid-toned, shortened surface is the low parameter estimate: The parameter is µ = 0.95, giving p ≈ √ 2. The points are data read from numerical solutions of the original Schnakenberg reaction-diffusion-convection equation (2) using NAG D03PCF.In the individual runs, the reactants are initially at steady state: (u, v) = (0, 0), with a small disturbance at x = 0.The boundaries are held at zero derivative: ux = 0, vx = 0.
the maximum estimate for large convection and strong diffusion on species v.However, for weak diffusion on v (δ < 0, ε < 1, log (ε) < 0), this maximum is not attained.
The general analysis remains viable: we can consider a different angle θ at the wave onset which gives other wave speeds.In particular, if θ = 0 we have a simple expression for the wave speed: γ + √ 2ε 1 − µ 2 .This value seems approximately correct for the case of weak diffusion on v (δ < 0, ε < 1, log (ε) < 0), with a smooth transition visible between this and the maximal estimate, as the diffusion on v varies.
Then we have a complete estimate for the wave speeds: with the right wave speed varying from the low estimate to high as the difference between the diffusion coefficients, ε, increases.

A. Increased reaction instability
We repeat the experiment for different values of µ.First we try µ = 0.85 (this gives α ≈ 0.14, β ≈ 0.84, p ≈ 1.43).As µ decreases the ODE becomes more unstable.For the PDE we expect little qualitative change to the system behaviour.This proves to be the case.The estimated wave speeds fit the numerical behaviour in the same way as the previous experiment.

B. Stable reaction, µ > 1
We extend the analysis to the case of α < 0. Our third experiment is for µ = 1.1, giving α = −1.05(and then β ≈ 1.09, p ≈ 1.42).Now we predict the special wave speed estimate, c .We expect We find the behaviour in our numerical experiment remains close to our estimate.In the combined case, significant difference between the two diffusion rates (ε 1) and α < 0, then the right hand side of the inequality could become positive again, and return to being a component of the wave speed estimate.This is effectively the requirements for the Turing instability.

VI. INVESTIGATING THE ANOMALY
We see that our Schnackenberg system behaves differently from our previously-analysed λ-ω system.We Biomath 2 (2013), 1303027, http://dx.doi.org/10.11145/j.biomath.2013.03.027 Figure 4.The θ behaviour (colour map) shows a broad propagation of pattern.However, the instability (r > 0, z-axis) does not propagate so widely.We see that the instability is only present for half of the underlying behaviour.The left wavefront is actually at the centre of the pattern, the line traced out from the initial disturbance.The parameters here are ε = 1.0, γ = 5.0, µ = 0.85.This is a numerical solution of the original Schnakenberg reaction-diffusion-convection equation ( 2) using NAG D03PCF.The reactants are initially at steady state: (u, v) = (0, 0), with a small disturbance at x = 0.The boundaries are held at zero derivative: ux = 0, vx = 0.
examine the behaviour of θ, and compare it to that of the radial parameter r.The θ diagram shows the system has an underlying propagation of pattern which is broad.We see that θ is behaving as it did in the λ-ω system.That is, it changes from being near-constant to increasing at a constant rate.This change occurs in a smooth, vshaped pattern, skewed to the right by the convection.However, the limit cycle (r ≈ 0) does not propagate so widely.That is, the behaviour of r is not directly related to that of θ.
By combining θ and r in Figure 4 we show that the instability is only present for half of the underlying behaviour.To the right, we see r instability linked to θ, initiating along the same line.This is as expected, and as the λ-ω system behaves.To the left, the onset of the wave is far from the θ initiation.The instability, or left wave "front", is occurring at the centre of the θ pattern.This is completely unexpected, and different to any behaviour we have seen before.
However, this different behaviour does correspond to the earlier fitting of the surfaces.Along this central line, θ is increasing continuously.This explains how the average value estimate for the wave speed continues to hold for the left wave speed.

VII. DISCUSSION
Reaction-diffusion systems produce many complicated patterns.A great deal of effort has been made to understand these systems.Here we examine a variation, a reaction-diffusion-convection system.We study the instability caused by a point disturbance from equilibrium.The disturbance expands in space linearly with time: it is consistently V-shaped.The advection carries the expanding pattern downstream.Within the disturbance, the chemicals oscillate.The oscillations vary in frequency and amplitude to the left and right of a line close to x = γ 2 t.In our previous work we analysed the behaviour of a simpler, symmetric system: the elementary λ-ω form [14].Here we look at a realistic chemical system.The reaction functions do not have the symmetry of the λ-ω system: we transform them so the linear component is symmetrical.This extension to a general system was suggested by Murray [15,]; here we demonstrate the application in a novel problem.From this, we find an effective way to predict the propagation of the instability.
We need to know the behaviour of θ to complete the prediction.Our analytical approach does not give this information; we turn to our numerics suggest some estimates.Our realistic model does not always behave like the simple λ-ω system.
The insight from our analysis sheds some light on the out-of-gamut pattern formation.That is, for µ > 1, the ODE is stable but the PDE in unstable.Our theory explains this: the pattern can form irrespective of the value of µ.This corresponds to the basic property of diffusion-driven (Turing) pattern.In our case, it is only the convection which appears to directly contribute to the propagation of the pattern.If convectivion is found in biological systems, then this system becomes directly relevant.
Convection-induced pattern formation complements Turing's idea of diffusion-induced pattern.Earlier work on this phenomenon showed some of the behaviour of the instability.Here we give some theoretical insight.Our results begin to explain the basis of the instability.

Figure
Figure2.Phase space for the limit-cycle reaction (3), given by numerical solution.The phase curve spirals out from the steady state to meet the limit cycle (broad loop).The Schnakenberg (a) system is highly elliptical, even at the start of the trajectory, by the steady state.The transformed Schnakenberg (b) shows a more circular behaviour initially, but the limit cycle shape remains more complicated.The parameter here is µ = 0.95.The solution was found using Matlab ode15s.

Biomath 2 (Figure 3 .
Figure 3. Numerical data compared to analytical estimates.The left wave speed is close to the low parameter estimate for all parameter values.In contrast, the right wave speed is far from the low parameter estimate, and fairly close to the high parameter estimate.While the behaviour corresponds to the analysis, the behaviour of θ is different to expectations.In the left wavefront speed, the darker surface is the low parameter estimate:γ/2 − √ 1 + ε 1 − µ 2 ;The lighter surface is the high parameter estimate:(1 − p)γ/2 − (1 + p + (1 − p)ε 1 − µ 2 .In the right wavefront speed, the lighter surface is the high parameter estimate:(1 + p)γ/2 + 1 − p + (1 + p)ε 1 − µ 2 .The similar dark surface is a less extremal estimate, with p cos φ = 1, giving the surface as γ + √ 2ε 1 − µ 2 ; the mid-toned, shortened surface is the low parameter estimate:γ/2 + √ 1 + ε 1 − µ 2 .The parameter is µ = 0.95, giving p ≈ √ 2.The points are data read from numerical solutions of the original Schnakenberg reaction-diffusion-convection equation (2) using NAG D03PCF.In the individual runs, the reactants are initially at steady state: (u, v) = (0, 0), with a small disturbance at x = 0.The boundaries are held at zero derivative: ux = 0, vx = 0.