Velocity-Amplitude Relationship in the Gray-Scott Autowave Model in Isolated Conditions.

Velocity and amplitude are two basic characteristics of any autowave, and their relationship reflects the internal regulation of the autowave system. This study proposes an approach to approximately estimate steady velocity-amplitude (VA) relation without deriving separate formulas for V and A. The approach presumes constructing an ansatz which represents the "petal" form of phase trajectory and contains V, A, and a free parameter (parameters). After substituting this ansatz, integration of model equations leads to a VA relation analytically. A free parameter (parameters) can be determined by comparing the analytical VA relation to the numerical data. As an example, we used the simplest autowave model possessing threshold, that is, the Gray-Scott model (cubic autocatalysis with linear inhibition) in isolated conditions with an immobilized precursor and a diffusible reactant. For all values of the inhibition rate constant allowing autowave solution (i.e., except approaching zero), the free parameter of ansatz belongs to a narrow interval has little effect on VA relation and can be regarded as fixed. Assumption of precursor immobilization does not influence the results qualitatively. This approach will be useful in investigations of more complex active media systems in biochemistry, combustion, and disease control.

Introduction

Active media can be generated from coupling positive feedbacks with mass transfer. They belong to the list of the most intriguing themes in chemical kinetics[1−4] and biology.[3−7] One of the most potent ways to theoretically investigate continuous active media systems is to use reaction–diffusion (RD) mathematical models of the least possible complexity. Such models give insights into the basic, qualitative properties of the system. Low-dimensional models can be either obtained by reducing large-scale detailed ones (usually comprising dozens of variables) or constructing phenomenologically from known/assumed general properties and/or (bio)chemical organization of the system. In this study, we focus on a special class of solutions of such equations, using the term “autowave” to refer to an active wave travelling in the form of self-sustaining pulse (with front, top and back), to distinguish it from a travelling front (backless).

Two basic characteristics of an autowave are velocity and amplitude. Generally, velocity and amplitude differently depend on parameters of reactions and diffusion. Thus, different systems should have different velocity–amplitude (VA) relations, and VA relation may be considered as a characteristic of the internal regulation of the system. This view poses the question of how VA relation can be derived from model equations.

The specific form of kinetic terms in RD equations (namely, nullclines’ shape, number and types of steady-states) largely determines properties of the system and suitable investigation approaches.[1−10] Nonlinearity of equations usually prevents their explicit analytical integration except for the simplest models which have travelling front solutions such as reduced Nagumo (or, purely cubic autocatalysis) and Kolmogorov–Petrovskii–Piskunov–Fisher (quadratic autocatalysis) models.[1,2,5−7,10] Interestingly, in the latter case, the explicit solution is known only for the special front velocity.[7,11] Even slightly more complex models are usually studied numerically or approximately. Approximate solutions are acceptable in mathematical biology and biophysics because models and their parameters are not exact, and qualitative properties of solutions are of greater interest.

Various analytical approaches have been developed for obtaining approximate solutions of RD equations and studying their properties: narrow reaction zone method (neglecting some reaction and diffusion terms in different space regions),[1] piecewise-linear approximation (replacing non-linear kinetic terms by some piecewise-linear function),[5−7] singular perturbation, or time/length scale separation analysis (investigating the limits of small parameter(s) tending to zero),[4,5,8,9] using ansatz (guessing the structure of a solution either in space or in a phase plane)[2,5,7,10,11] and their combinations. In most cases, these approaches were applied to study travelling fronts and autowaves with saturated amplitude. So far, general methods to estimate the VA relation for autowaves do not exist, to the best of our knowledge.

In this article, we propose the new, to the best of our knowledge, analytical–numerical method for estimating the VA relation in autowave systems. The method is based on approximating the phase trajectory of the system by an ansatz—single smooth function rendering solutions in all space regions and containing V, A, and a free parameter (parameters). After substituting this ansatz, integration of model equations leads to a VA relation analytically. A free parameter (parameters) of the ansatz can be obtained by comparing this VA relation and direct numerical solution data. As an example, we take the simplest autowave model with threshold properties: a two-component system with precursor (v) conversion to the autocatalyst (u), obeying cubic kinetics, followed by a linear inhibition of the autocatalyst, see Scheme . Sometimes it is referred to as the Gray–Scott (GS) model. We suppose that the system is isolated (i.e., the influx of components is absent).

Scheme 1

Kinetic Scheme of the two-Component System with Cubic Autocatalysis and Linear Inhibition (GS Model)

Various complex systems can be studied after reducing their kinetic scheme to the GS model.[2] Examples include glycolysis,[12] catalytic chemical reactors and cold flame propagation,[1,9,13−17] fungal mycelia growth,[18] and disease spreading.[19] However, available studies in the spatial-distributed conditions are almost always limited to open conditions,[9,17−20] the dependence of autowave velocity on parameters in isolated conditions is rarely studied,[13,16] and the VA relationship in the GS model in isolated conditions has never been reported, to the best of our knowledge. Here, we study this relation by the proposed analytical-numerical approach which, as we suggest, can be further adapted for more complex autowave models.

Description of the Model and the Analytical Approach

Model Assumptions

We consider the GS model (Scheme ) under the following assumptions:

Geometry and boundaries: The system is one-dimensional and isolated (there is no volume influx of components), and boundary influence is absent or negligible.

Reaction kinetics: All reaction kinetics are subjected to the law of mass action. Inhibition rate of the autocatalyst is nonzero, so the shape of its concentration profile is pulse (see Figure A, thick line).

Figure 2

Numerical solutions and analytical orbit approximations in the GS autowave model. (A) v(x) and u(x) profiles found by solving eq at h = 0.04. Steady autowave velocity and amplitude are c = 0,70 and m = 0,725. From eq , λ1(c,h) ≈ 0,053 and λ2(c,h) ≈ −0,75, and from eq , I(c,m,h) ≈ 2.90 and n2(c,m,h) ≈ 1.31. (B) Same numerical data in the u–u′ plane (black solid line), their best fit by eqs and 9 with λ1 ≈ 0,054 and n1 ≈ 5 at p > 0 and λ2 ≈ −0,72 and n2 ≈ 1.27 at p < 0 (gray solid lines), effect of ±5% variation of n2 (dashed half-orbits), and tangents given by eq with the same λ1,2 (dash-dotted straight lines). (C) Numerical solutions of eq at three h-values (solid lines). Analytical half-orbits (p < 0) given by eqs and 9 with λ2(c,h) from eq and n2(c,m,h) from eq at the same h-values (dashed lines).

Initial conditions: Initially, the precursor is uniformly distributed throughout the space, and the autocatalyst is absent everywhere except the “seed” quantity localized near one of the boundaries (for Figure A, near the left boundary).

Steady-state regime: The system has enough time to transit to the steady travelling wave regime, and transient process is not considered in this study.

Diffusion: Precursor diffusion may be neglected as qualitatively unimportant process provided the ratio of diffusion coefficients of precursor and autocatalyst is not large. Actually, the main effect of precursor diffusion is some decrease of its concentration compared to the initial finite value in front of the autowave front (in contrast to the autocatalyst diffusion which enables medium activation and wave propagation). The effect of this simplification will be addressed at the end of this article (see the Discussion and Figure A,B). This assumption is a fortiori valid for systems with an immobilized precursor.

Figure 4

Numerical and analytical VA relations in the GS autowave model. (A) Dependences of autowave amplitude and velocity in the numerical solution of eq on the inhibition rate constant (solid lines). Results with added v term to the equation for v are plotted by dashed lines. Dotted line shows c(h) in the reduced Nagumo approximation, eq . (B) Same numerical data in the c–m coordinates (bold solid and dashed lines). Thin gray lines show the analytical c–m approximation, eq , with parameters: dotted line, (h,I) = (0.01,3.7); solid line, (h,I) = (0.04,2.9); dash-dotted line, (h,I) = (0.07,2.6). The area formed by the first and the last of these (h,I) pairs is depicted by rectangle in Figure B.

Governing Equations

The standard form of RD equations corresponding to the Scheme in isolated spatially-distributed conditions and at neglected diffusion of v, iswith boundary conditionsand initial conditionsHere, u and v are concentrations of the autocatalyst and its precursor, k and h are the activation and inhibition rate constants, D is the diffusion coefficient of the autocatalyst, v0 is the initial uniform precursor concentration, and u0(x) is the initial autocatalyst concentration, which is supposed to be nonzero only in the small but finite spatial domain. If both u and v diffuse, the equation for v in eq should be appended by the D·v term.

The nondimensional form of eq can be obtained defining t̃ = t/kv02, ũ = u/v0, ṽ = v/v0, , and h̃ = h/kv02. Omitting tildes gives

Boundary conditions remain the same, and δ0(x) ≈ 1 in the finite spatial domain and zero elsewhere. Here, F(u,v) = v·u2 – h·u. Note that F′(u = 0) = −h < 0, which gives threshold properties to the system.[8,21] Typical solution of these equations is shown in Figure A. Finally, to study solutions travelling with a steady velocity, c, that is, stationary in the moving frame x̂ = x̃ – c·t̃ so that (u,v)(x̃,t̃) = (u,v)(x̂), eq can be rewritten as

In this eigenvalue problem for c, prime (′) denotes the derivative with respect to x̂, and below x̂ will be referred to as x.

VA Relation from the Homoclinic Orbit Approximation

In this article, we suggest the following approach to estimate the VA relation for autowave systems that are qualitatively similar to the GS model. Let us denote and m = max(u) in the solution of eq . In the (v,u,p) phase space, this solution is a heteroclinic orbit connecting steady states (1,0,0) and (v2,0,0), corresponding to x → +∞ and x → −∞, respectively. Projection of the orbit to the (u,p) plane is a homoclinic orbit starting from (0,0) and returning to (0,0), see Figure B, for example. The homoclinic orbit intersects the 0u axis in the point (m,0) with the slope . If we construct the formula (ansatz) approximating the shape of the orbit and containing c and m, then after integration, we will obtain an approximate VA relation.

(a) At x → ±∞, we have: u → 0 and vu2 ≪ hu. Thus, u approximately obeys the homogeneous equation not containing v

The characteristic equation and its roots are

Solution of eq is

Thus, the orbit at the origin has tangents (in the (u,p) plane, see straight lines in Figure B)

(b) At 0 ≤ u ≤ m, we suggest that the ansatz approximating half-orbits at p > 0 and p < 0 has the formso that these halves smoothly join at the point (m,0). To do this, we need functions g(u) with properties

The last condition guaranties smoothness of the junction. In this study, we use the specific form of g(u) that fulfills all these conditions

At n = 2, this equation defines an ellipse arc (Figure A). We assume that n > 1, and that n may depend on h. With the above g(u) form, the equation for v in eq can be integratedwhere v = v1 at u = m. The integral I can be takenHere, , n ≡ n2, Β and Γ are beta and gamma functions, I1 = ∞, and I2 = 1. Equation uniquely binds n and I (Figure B, solid line). In practice, it is convenient (but not necessary) to use a simpler approximate formula. All I values in this study are well below 10 (see Figure B, gray line), and for I ∈ [1,10], the following relation takes place

Figure 1

(A) g(u) function (eq ) at different n values. (B) I found numerically from eq (black solid line) and approximated by eq with a = 0.15 (gray dashed line).

Figure 3

Parameters of the GS autowave model as functions of the inhibition rate constant obtained by different methods. (A) Parameters of the bottom half-orbit λ2 (left axis, gray) and n2 (right axis, black). λ2 and n2 (closed squares and closed circles) from the best fittings of the numerically obtained orbits at simultaneously varying λ2, n2, and m. n2 (open circles) from varying of this parameter solely (while λ2, eq , and m were taken from numerical results). λ2 and n2 (gray and black lines) calculated by eqs and 18 from the numerically obtained c and m, which are plotted in Figure . (B) I(h) and n2(h) (gray and black lines) found from numerically obtained c and m (Figure ) using eq . For comparison, n2 found by fitting of orbits is also shown (closed circles). Rectangle shows h and I bounds used in Figure B to plot the analytical VA relation.

(Figure B, dashed line). The value a ≈ 0.15 was obtained by numerical minimization of the integral of the absolute difference between I given by eqs and 12 for n ∈ [n*,2], where n* ≈ 1.093 is the numerical solution of the equation I(eq 11) = 10.

(c) Another relation between v1 and m can be found by considering the equation for u in eq at u = m. At this point, u′ = p = 0, v = v1, and g = 0. To find the value u″ = p′p = ∞·0, we should calculate

Substituting u = m and g(m) = 0 gives p′p = 0 – λ2m·g(m)2–, which is 0 if n < 2. Thus, eq becomes the simple equality of the inhibition and production rates

(d) The c–m relation can be found by equaling the expressions for ln v1 from u- and v-kinetics, namely eqs and 10

This relation can be rewritten as

From eqs and 12, we can find I and n2 as functions of autowave velocity and amplitude

Reduced Nagumo Approximation

Another way to analyze eq is to reduce them to a simpler model which is unique in that it has exact analytical solution. This model is known as the cold (isothermal) flame model,[1,13] purely cubic autocatalysis model,[2,10] or reduced Nagumo model.[5,6] In the region of steepest gradients, the rate of v → u reaction prevails over those of u inhibition and diffusion, thus v + u ≈ 1. More accurately, this first integral can be obtained by summation equations for u and v in the limit D/D → 1 and h → 0. Setting v ≈ 1 – u giveswhere

Equation describes the travelling front connecting two stable steady states, 0 and u3. This equation can be solved explicitly by guessing that the phase trajectory is u′ = p(u) = A·u·(u – u3).[1,2,5,7,10] Substitution of this guess into eq gives that , and front velocity equals to

Results

Overall View of Numerical Solution

The representative numerical solution of eq (for h = 0.04) is shown in Figure A. The autowave travels from left to right with steady velocity and shape. The v(x) profile has a switch-off type (thin curve), and the u(x) profile has a pulse type (thick curve). The same data in the u–p plane are shown in Figure B (black solid line), and solutions at several h-values are also shown in Figure C (solid lines). All solutions plot nearly homoclinic orbits starting from (0,0) to the bottom quarter-plane, intersecting the 0u axis from below and tending to (0,0) from the top quarter-plane. Ingoing of orbits to (0,0) ceases at low h (Figure C) because of finiteness of spatial interval of numerical integration, but this region of wave tail does not influence wave behavior, and increase of spatial interval should make all orbits definitely homoclinic (not shown).

Comparison of Numerical Solution and Analytical Orbit Approximation

Surprisingly, an ansatz constructed in eqs and 9 gives very good approximation of numerical orbits for every h when parameters λ, n, and m are varied just to get the best fit. Gray solid lines in Figure B show fitting example for h = 0.04, and dash-dotted lines show tangents (eq ). Parameters λ2 and n2 of these fittings at different h values are plotted in Figure A (closed squares and closed circles). Obtained λ2 and m are close to those calculated from numerical solutions (see gray squares vs gray lines for λ2). That is why similar fittings were obtained when only n2 was varied but λ2 and m were fixed to values obtained from numerical solutions (Figure A, open vs closed circles). Thus, the constructed ansatz represents the shape of orbit well and may be used to study general properties of the system.

Sensitivity of the Orbit Shape to the Parameter n2

Equations and 9 mean that the orbit trajectory at x → ±∞ and u → m is mostly linked to velocity (through λ) and amplitude, respectively. To check that n influence is localized only in the intermediate zone, two dashed half-orbits in Figure B show the effect of ±5% change of the n2 value. Indeed, the effect of variation of n2 is maximal between points (0,0) and (m,0), in the region of the steepest u(x) gradient. However, there is no any “hypersensitivity”, as 10% change of n2 leads to ∼11% change of |p|max.

Ansatz Shape Parameter n2 can be Estimated from the Numerical Velocity and Amplitude

The next aim is to obtain the ansatz shape parameter related to autowave velocity and amplitude using model equations (rather than by orbit fitting). Numerical values of c and m in the entire range of h values are plotted in Figure A (solid lines). Expectedly, the steady autowave solution only exists at h below some threshold value, h*. We have found that h* ≈ 0.07. At h > h*, the wave does not survive irrespective of the triggering stimulus, δ0(x), and c = m = 0 are plotted. The same data in the c–m coordinates are plotted in Figure B (bold solid line), considering h as a parameter. Because of threshold, the c–m line does not reach the origin. Assuming that the derived analytical VA relation (eq ) may be applied to these data, I(h) and n2(h) were found using eq (Figure B, gray and black solid lines). The principal difference between these n2(h) (black line) and n2 found above by simply fitting the u–u′ half-orbits (circles) is that we derived eqs –18 taking into account model equations for both u and v, and that we used the computed autowave velocity in eq . Half-orbits with parameters λ2 and n2 calculated from numerically obtained c and m (eqs and 18) are plotted in Figure C (dashed lines). They do not match numerical orbits precisely but at x → +∞ and at u → m they do.

Obtained I(h) and n2(h) dependences are sloping quite mildly except for the region h → 0 (Figure B). Because at h = 0, the u(x) profile in solutions of model equations has a shape of travelling front but not pulse, below we only consider nonzero h values.

Analytical VA Relation Weakly Depends on Particular (h,I(h)) Values

Figure B shows that at almost all h (except h → 0), the values of I(h) and n2(h) are bounded in the narrow ranges. For example, 2.6 ≤ I ≤ 3.7 and 1.24 ≤ n2 ≤ 1.34 at 0.01 ≤ h ≤ 0.07 (shaded rectangle). Substitution of these h and I bounds into eq gives the narrow corridor bounded by the dotted and dash-dotted gray lines in Figure B, which contains almost all numerical c–m results (bold solid line). Thus, any line passing through this corridor is an appropriate approximation of the c–m relation in the selected range of h values. For example, solid gray line in Figure B is plotted for intermediate values (h,I(h)) = (0.04,2.9). The deviation is small even for h → 0 (m → 1 in Figure B).

Discussion

VA relation is an important but underestimated characteristic of any autowave system. Usually, velocity is the only autowave property under consideration: first, nonzero velocity shows the very possibility of signal propagation, and its value gives the time required for a signal to cover the distance; second, the biological system response is usually regarded as a Yes/No (All/Nothing) type, so the precise value of autowave amplitude looks somewhat unimportant. However, autowave amplitude can be of major interest in various situations, such as the following:

The main outcome of a system is determined by autowave amplitude. For example, in blood clotting, the thrombin autowave leaves behind the polymerized fibrin mesh.[22,23] The fate of clot (to reside competing with blood flow or to embolize partly or wholly) depends on its physical, that is, rheological, properties such as storage and loss shear moduli. These macroscale properties depend on microscale characteristics of fibrin mesh (fibers’ diameters, lengths, density, branching density, and cross-linking), which are determined mainly by the local transient thrombin concentration.[24−26]

Autowave amplitude has great value in the economic, humanistic, and so forth, sense. For example, it designates the number (or fraction) of infected persons during spatial spreading of a disease.[6,19]

The amplitude is easier to measure then the velocity, or the velocity is too low/too large to measure it correctly. In this case, the VA relation (like eq for the GS model) can be used to estimate the velocity from the measurements of amplitude.

In this study, we propose an approach to investigate VA relation in autowave systems in steady travelling regime. The key distinction of this approach is that it does not regard derivation of separate formulas for velocity and amplitude as functions of model parameters. Instead, it presumes the following steps:

Approximating the orbit shape in the phase space by a “petal” or “arc” nonlinear function (ansatz), which contains wave velocity and amplitude as well as a free parameter of orbit curvature.

Substituting the ansatz into model equations and integration of equations all along the wave front to the wave maximum; after some algebra, this gives an analytical VA relation.

Calculating dependence of the curvature parameter of the ansatz on model parameters from the derived analytical and numerical VA relations.

Checking whether the curvature parameter of ansatz has little effect on VA relation and may be regarded as fixed.

In our study, the free parameter, n2, and its function, I (eq ), turned out to be bounded within the narrow ranges for all values of the model parameter h, allowing travelling wave solution (i.e., except for h → 0), see rectangle in Figure B. Any (h,I(h)) pair from this range can be used in the analytical VA formula (eq ) to well match the numerical results for all h (Figure B). This agreement in VA relation did not require perfect matching between analytical and numerical orbits (see Figure C), suggesting that the orbit shape is the “finer” property of the system compared to the more “robust” VA relation property.

Moreover, the orbit shape may be approximated in numerous ways, and the choice of the function g(u) looks arbitrary. Indeed, g(u) in eq may be multiplied by, for example, 1 – αu2 or exp(−αu2), still fulfilling the required conditions (eq ) but introducing an additional free parameter, α. For the model system considered in this study, this complication looks unnecessary. Another advantage of the chosen g(u) form (eq ) is the integrability of I (eq ), although this is not strictly necessary for the method, as eqs and 17 contain I but not n2. The question of the preferred choice of g(u) functions at different kinetic schemes will be addressed in a separate study.

Ansatz construction is a powerful but unformalized approach to obtain approximate and exact solutions of nonlinear differential equations.[7] The idea is to “guess” the structure of a solution basing on the knowledge of its general properties, and to test this guess by substituting it into the model equations. For example, the generalized Kolmogorov–Petrovskii–Piskunov–Fisher (KPPF) equation, that is eq with F = u·(1 – u), can be solved by constructing an ansatzu(x) inspired by the main term of the approximate solution of the “classical” KPPF equation having q = 1.[7] Another example is suggestion of the polynomial form for the heteroclinic phase trajectory, p(u), if the reaction term in eq is polynomial.[7] For cubic kinetics, this means that the phase trajectory can be the second order parabola, which is often used as a “guess”, allowing an exact solution.[5] However, for autocatalytic systems with two variables (Scheme ), this approach can only be used in case of approximately equal diffusion coefficients and zero inhibition rate, allowing to reduce a model to one equation (like eq ).[1,2,10,13] Thus, all that has been said above applies to the models which are reducible to one equation and have front-like solutions. So far, there has not been constructed an ansatz for a model with two variables and nonzero inhibition rate of the autocatalyst.

Cubic autocatalysis with linear inhibition is a well-known model system arising in numerous chemical and biological tasks, from flame propagation to glycolysis. Even in the homogeneous case, this model was shown to have “exotic” (for chemistry) properties such as bistability and time-periodic solutions, provided there exists a continuous influx of the first component, that is, the conditions are open.[12,14,15] After the last articles, this model is sometimes called the GS model. The spatially distributed GS model was studied in several articles,[9,17,18,20] again in open conditions. In short, Petrov et al.[17] studied regimes of autowave propagation, reflection from boundaries, and from each other in the presence of influx of both components; Davidson et al.[18] studied fungal mycelia growth at the continuous influx of v (replenishment of the substrate); Hale et al.[20] derived the formula for the orbit at the large influx rate of v, equal to (v0 – v)hD/D (in notations of eq ); Muratov and Osipov[9] used the singular perturbation technique to estimate autowave amplitude and velocity at the nonzero influx/removal rate, k0; and their solutions diverge at k0 → 0. Equations similar to the GS model emerged in the susceptible-infective model with the nonlinear infection rate,[19] with the only difference that the precursor influx is replaced by the reproduction rate of susceptibles. However, in isolated conditions, the GS model was considered only in two available studies.[13,16]

Novozhilov and Posvyanski theoretically studied the cold flame propagation experimentally observed earlier by Voronkov and Semenov in the 0.03% CS2 burning in the excess of oxygen, outside the ignition peninsula.[1,13] Considering the branched chain reaction scheme, they derived RD equations similar to the eq (but with both components diffusible) and solved them numerically. Merkin and Needham[16] studied the same system analytically and numerically. In both studies, solutions existed only for h < h**, and the c(h) dependences were similar to that obtained in this study (Figure A, gray solid line). To make the direct comparison, the v term was added into the equation for v in eq . Results of numerical integration of the “extended” system are also shown in Figure A (dashed lines). Solution exists for h < h** ≈ 0.045, in agreement with h** = 0.04611[13] and 0.0465.[16] At the zero inhibition rate, velocity c(h = 0) ≈ 0.7 also agrees with the published results[13,16] and with the reduced Nagumo approximation, (eq , dotted line in Figure A).

Comparison of solution with and without diffusion of v (Figure A,B, dashed vs solid curves) shows that diffusion of v, at least in the range 0 ≤ D/D ≤ 1, affects properties of the system quantitatively but not qualitatively: it reduces autowave velocity and amplitude as well as h* but does not disrupt V(h), A(h), and VA relations. Novozhilov and Posvyanski obtained similar V(h) curves when varied the ratio D/D in the range 0.5–2.[13] Qualitative unimportance of diffusion of v (at not very high ratio D/D) can be well understood as it acts to smooth the v(x) profile: to decrease v at the wave front (where v ≈ 1) and to increase v in the wave tail (where v ≈ 0). Some decrease of v compared to v ≈ 1 at the wave front can indeed decrease the velocity, but the wave tail does not influence the system behavior. On the contrary, diffusion of u is a qualitatively important process as it brings u to the region where u ≈ 0, and the medium is activated in front of the front (see Figure A, x ≳ 85). Relative (qualitative) unimportance of v diffusion is the main reason we neglected it in this study. Another reason is that it makes equation for v ordinary and thus easier to integrate. Similar simplifications are used in other systems with low-mobile precursors, such as disease spreading (for susceptibles),[6,19] predator-prey (for pray),[6] resource-consumer (for resource),[7,27] and so on.

Numerical Methods

To integrate eq numerically, we used the Störmer–Encke method for space discretization (uniform mesh with N = 401 nodes, zero Neumann boundary conditions at x = 0 and x = 100) and the embedded Runge–Kutta–Fehlberg method of order 2(3) with automatic step size control for integration in time (mixed local error estimation with max norm, tol = 0.01, fac = 0.8, facmax = 5).[28] Increasing N and decreasing tol had negligible influence on c and m.