FRONT PROPAGATION PROBLEMS WITH SUB-DIFFUSION

Front propagation is considered for two kinds of sub-diffusion – reaction systems: (i) systems with sub-diffusion limited reaction rate governed by models with fractional time derivatives; (ii) systems with activation limited reaction rate governed by integro-differential equations with two time variables. It is shown that in the former case the front is described by a travelling wave solution, while in the latter the velocity of the front decreases with time.


1.
Introduction. Anomalously slow diffusion, alias sub-diffusion, is a dispersion of particles whose mean square displacement is given by a power law r 2 (t) ∼ t γ , 0 < γ < 1. The limit of normal diffusion is obtained with γ = 1.
Sub-diffusion often characterises the motion of particles within living cells due to the complicated structure of the medium. One example is the diffusion of enzymes within a cytoplasm that is a moderately crowded environment, whence a typical anomaly exponent γ is expected to be close to unity [24]. The motion of transcription factors within nuclei is severely hindered by the closely packed DNA strands, and the anomaly is more pronounced with the expected value of γ about moiety [23]. The lateral diffusion of molecules along cell membrane is combined from a rather mobile fraction ( e.g. ligands ), whose anomaly exponent is moderate, and an almost immobile fraction ( e.g. receptors ) with γ → 0 [7,15]. In all of these examples and many others the particles travel through the cell and its organelles in order to complete a destined reaction.
The reactions can be diffusion limited, when the reaction rate is determined mainly by the availability of reactants or catalysts, or activation limited if the energy supply is precluded. In the former case due to the anomalously slow diffusion the reaction rate is generally diminished. This phenomenon was demonstrated in one experimental work with yeast cells [20] and a few theoretical studies of sub-diffusion with different reactions: coagulation and annihilation [26], front propagation [8], reaction with a singular disorder [5], reaction with bimolecular enzyme kinetics [6] and reaction with trapping in exciton dynamics [27]. There are exceptions, however, and an enhanced reaction rate was obtained in diffusion limited reaction due to anomaly [16]. This is to show that the diffusion anomaly might result in nonintuitive changes of the reaction rates. In the case of activation energy limitation the reaction rate might stay unchanged.
Due to the wide variety of anomaly origins in diffusion processes within living cells, no universal model of sub-diffusion can be formulated, and the mathematical modelling of sub-diffusive reaction processes is far from being complete. By a recent study, diffusion anomaly on a cell membrane occurs by three different mechanisms [14]. Sub-diffusion has been traditionally modelled by fractional derivatives [11]. In fact, there are many types of fractional derivatives and moreover, it has been shown that certain reaction -diffusion processes with anomaly can only be modelled by more complicated integro-differential operators [19,17].
One of the pertinent reaction -sub-diffusion problems is that of front propagation. Some reactions in living cells are front-like. Among these are signals sent between the cell organelles [4], fronts travelling during tumour invasion into healthy tissues [9] and fronts produced during activation of actin filaments [21].
Sub-diffusive fronts have received scarce attention due to the aforementioned modelling difficulties. Inspired by the numerous applications in cell biology, in the presented work a basic theory is set for front propagation subject to sub-diffusion with an emphasis on the differences due to the choice of the memory operator. In section 2.1 an expression for the propagation velocity is obtained for a fractional derivative as the memory operator and a piecewise linear kinetics. In 2.2 domain walls are analysed for Allen-Cahn kinetics and a small external field. In section 3 a similar analysis is performed for a more complicated memory operator, involving the concept of aging.
2. Models based on fractional derivatives. Hereinafter two fractional subdiffusive models are considered, touching on different aspects of front propagation. The proposed equations describe a sub-class of reactions and because of the diversity of anomaly origins in nature are not meant to model a specific system.
First, a simple fractional model is analysed, amenable to an exact solution due to a piecewise linear kinetic function. This basic, single dimensional model depicts how the presence of anomaly complicates the dependence of the front speed on the system parameters. Next, a fractional model with a cubic non-linearity and a small external field is investigated. Here changes of front profiles and velocities are studied.
A normal reaction -diffusion equation reads wherein u is the species density, D -a constant diffusion coefficient and f -a kinetic function. One class of sub-diffusive problems is obtained via a modification of (1) by a fractional derivative of order γ where the fractional operator is defined as [11,3] wherein a kinetic function was chosen so that (4a) is amenable to an analytical solution with k > 0 being a kinetic constant, H denoting the Heaviside function and a -an arbitrary discontinuity point ( see figure 1 ). This choice of the kinetic function allows for an exact solution in the form of a travelling front, i.e. satisfying (4c) The section will focus on the propagation of the front and not its formation, thus enabling the assumption that the front has existed indefinitely long ago and using an infinite lower bound in the memory operator (4b). The velocity c is so far unknown and might be of any sign. Due to the translational invariance of such a front it is possible to assume and then using the boundary and monotonicity conditions (4c) replace H(u − a) by H(x).
2.1.1. Front shape. Fourier transform of (4a) with respect to ξ yields with q being the transform variable and transformed functions denoted by a hat. The transform is taken in the sense of distributions, and δ is the Dirac delta function.
With the aid of lim equation (6) is inversed as The integral can be readily computed for the part of the wave that is 'ahead of the front', i.e., for ξ > 0 in the case when the wave propagates to the right (c > 0), and for ξ < 0 in the case when the wave propagates to the left (c < 0). To evaluate the integral, the residue theorem is used with the contours constructed from the integration path along the real axis and a semi-arc of radius R in the upper half plane for c > 0, ξ > 0 and the lower half plane for c < 0, ξ < 0 ( see figure 2 ). As the limits ǫ → 0 and R → ∞ are approached, the contour integral equals the desired one, since the arc integral vanishes for both upper ( c > 0 ) and lower ( c < 0 ) contours, as seen by substituting q = R e iθ to yield Thus It is assumed that the complex function (−icq) γ has a branch cut along the positive imaginary axis for the case c < 0, ξ < 0, and along the negative imaginary axis for c > 0, ξ > 0.
Substituting q = sgn c iQ, Q > 0, the integrand possesses a pole Q ⋆ on this branch for c ≷ 0 unique over R by a graphical solution method ( see appendix B for the proof of uniqueness over C ). For c > 0 there is one more pole q = iǫ within the integration contour. Therefore, bearing in mind that for c < 0 the winding number of the contour is -1, For c > 0 the winding number of the contour equalling unity. It is seen, in particular, that ahead of the front the solution approaches its limiting behaviour at infinity exponentially fast. The behaviour of the solution behind the front is qualitatively different from that ahead of the front, in particular, it does not behave exponentially. For sgn c = sgn ξ it is convenient to write the solution (8) in the real form as sgn c |c| γ q γ sinγ cos(qξ) (D q 2 + k + |c| γ q γ cosγ) 2 + |c| 2γ q 2γ sin 2γ + D q 2 + k + |c| γ q γ cosγ sin(qξ) The asymptotics is obtained via the integral form (13) and the decay is algebraic ( see appendix C for a detailed derivation ) The difference in the functional form of the tails corresponds to the system's memory, as induced by sub-diffusion.

Propagation velocity.
Imposing u(0) = a it is seen immediately that (12) is continuous at c = 0 for all 0 < γ 1. Combining with (11) and solving for c, Thus, c is an odd function of a − 1/2 that is equal to zero at a = 1/2 and monotonically increases with a. As a goes to 1, c goes to infinity, In the case of normal diffusion (15) reduces to Consider, e.g., (18a). If a > a cr ≡ k + 2 2(k + 1) , so that c is extremely large, i.e., propagation of the front in the case of strongly anomalous diffusion is much faster than that in the case of normal diffusion. If, on the other hand, 1/2 < a < a cr , the propagation speed is extremely small. This non-intuitive threshold behaviour can be explained as follows. Normally, condition (5) implies that if c = 0 ( a = 1 2 ), reaction and diffusion terms balance each other, and the front is stationary. If c = 0, the front moves so that u = a on the line x/t = c in the (x, t) space. When the state u = 1 ( a > 1 2 ) dominates, a rightward motion ( c > 0, x > 0 ) yields the desired balance, and vice versa when u = 0 dominates. With a > 1 2 and a slight anomaly the balance is reached despite the slown diffusion, whereas for a sufficiently large the front must speed up to "catch up" with the reaction.

2.2.
Fractional model with a cubic non-linearity. In order to show that the front properties found in the previous subsection are typical for sub-diffusionreaction system, problem (2b) is now considred with the fractional operator defined by (2c) and a smooth kinetic function f (u) = u − u 3 + µ corresponding to the Allen-Cahn equation [1]. In the case of normal diffusion that equation describes the kinetics of a phase transition. The variable u is now the deviation of the density from a certain characteristic density rather than density itself, hence it can take both positive and negative values. Parameter µ, which breaks the symmetry u → −u, plays the same role as 2a−1 in the piecewise linear model considered in the previous subsection. Similar models with different kinetics were studied in the past [10].
A solution of the equatioñ is sought in the form of a front with velocity c, The boundary conditions ( stated below ) are to be chosen so that for µ > 0 the domain wall moves to the left, i.e. c > 0. Then if µ < 0, the transformation µ → −µ, u → −u, x → −x restores the considered problem. Hence the front solution is governed by the equation The boundary conditions are taken as where u + (µ) and u − (µ) are the maximal and minimal roots of the equation ( it is assumed that |µ| < 2 √ 3/9, therefore equation (22) has three real roots; note a difference in sign in the definition of c in comparison with the previous section ). In the case µ = 0 equation (19) has a stationary solution u 0 (ξ), c = 0, which satisfies solved as This is similar to the case a = 1/2 for the model considered in section 2.1. Let us consider now the case 0 < µ ≪ 1 and construct the solution in the form At the leading order, the correction u 1 satisfies For (23c) to have a solution, by Fredholm alternative with the quantity b(γ) given by The case µ < 0 can be considered in a similar way; one finds that c(−µ) = −c(µ). The obtained relation for c(µ) is similar to the relation (15) found in the previous subsection.
The effect of memory is seen when the asymptotic form of the solution tails is evaluated. At the left tail ( ξ → −∞ ) the fractional term on the left-hand side of (23d) yields and the asymptotics for u 1 is with C being a constant that cannot be determined by this approach. At the right tail ( ξ → ∞ ) the substitution ζ = yξ is used in the evaluation of the fractional derivative, giving The obtained asymptotics are similar to those found in the previous section. The difference of the asymptotics for ξ → −∞ and ξ → ∞ can be explained qualitatively in the following way. As the front propagates to the left, the left tail is not "aware" of its existence and exhibits a normal exponential decay due to the hyperbolic tangent solution, whereas at the right tail the memory of the already passed front causes the power law asymptotics typical of the kernels of fractional derivatives.
3. Domain walls for sub-diffusion with aging.
3.1. Models of sub-diffusion based on continuous time random walk approach. In the present section another class of sub-diffusion models is addressed, where the rate of chemical reaction is activation limited and not influenced directly by the sub-diffusion. Due to the presence of memory the reaction and diffusion processes are intertwined in a way that no longer allows to write a balance equation with the respective terms uncoupled. The formalism of continuous time random walk is utilised to build a positive definite evolution operator that accounts for sub-diffusion with reactions [19,17,22]. Let us describe first the essence of that approach in the absence of a reaction. Assume that the system contains N kinds of molecules performing random jumps starting from a certain time instant t = 0. A molecule of the kind i, i = 1, 2, . . . , N , that arrives to the point x ′ at the time instant t ′ performs the next jump to the point Let us denote by n i (x, t, t ′ ) the density of molecules of the kind i at the point x and time instant t, which have arrived to that point at a time instant t ′ , 0 < t ′ < t.
The time evolution of that quantity is given by hence Using (29), it is possible to rewrite (28) in the form where Equation (30) is valid for t > 0. To account for the initial condition the definition of functions n i (x, t, t ′ ) is formally extended into the whole region −∞ < t < ∞, so that n i (x, t, t ′ ) = 0 everywhere outside the subset 0 t ′ t. At the instant t = 0 particles of densities ρ 0i arrive in the system, so that n i (x, 0, t ′ ) = ρ 0i δ(t ′ ). At finite t the distribution n i (x, t, t ′ ) consists of two parts: a smooth part at t ′ > 0 comprising particles that performed jumps, and a singular part of the type ρ si (x, t)δ(t ′ ), where ρ si (x, t) is the density of particles that never performed jumps. The total probability density is wherein the limit 0 − of the integral means that the contribution of the δ-like part of the distribution is taken into account. At the same time, the general master equation reads Comparing (28) and (30), one comes to the conclusion that equation (33) can be written also as Later on, w(τ ) is taken as giving A "heavy tail" in the waiting time probability function provides a sub-diffusive behaviour [11]. Let us include the chemical reaction. Now, the time evolution of n i (x, t, t ′ ) is affected by two factors: total decrease due to jumps and chemical composition change due to the reaction. The following modification of equation (30) is postulated [13]: with the probability density for walkers of all ages being (37) It is convenient to introduce the variable τ = t − t ′ ( the particle "age" ) instead of the arrival time t ′ [22] and define Hence the model is formulated as For the sake of simplicity, later on the limit t is indicated instead of t + in the integrals. It is meant however that a δ-like part of the distribution stemming from the initial condition is taken into account. The model desribed above is based on the assumption that (i) each kind of molecules is produced from all kinds of molecules with the rates M ij that do not depend on the ages of molecules but only on the total densities of components; (ii) the age of the molecule does not change in the course of reaction. The latter assumption needs some commenting.
When a diffusive process is described by a continuous time random walk, a particle's age is the time elapsing between its last step and the present time. If the process also involves chemical reactions, the particle might change identity before its next step occurs. Then two possibilities ensue: the definition holds and the particle keeps its age in spite of the change of its identity, or the definition is adjusted so that any newborn particle is ascribed the age of zero. Both approaches were studied in the literature, and there are arguments in favour [13] and against [25] each one of them. On one hand, it is more intuitive to assign a zero age to a new particle and quite difficult to track all the particles of the same age, when mutual conversions might have taken place more than once. On the other hand, the "rejuvenation" of particles in the course of the chemical reaction, even near the state of a chemical equilibrium, may significantly enhance the diffusion and make it insensitive to the behaviour of the tail of the waiting time distribution w(t), because "young" particles jump more often, and their rate of jumping does not depend on the tail asymptotics [2]. In the present paper, the front-type solutions for a sub-diffusive model with uniform aging are considered.

Spatially uniform solution.
Prior to the discussion of the front propagation, let us consider a spatially uniform solution governed by the following system of equations and boundary conditions: where the integrals include contributions of both the smooth distribution with τ < t and the singular distribution at τ = t. The general solution of (40a) is and by (40b) Of course, in the spatially uniform case the diffusion is absent, therefore the age distribution of molecules is irrelevant. Indeed, substituting relation (41b) into the expression for ρ, one obtains a closed equation for density, Thus, for sufficiently large t ρ(t) tends to a constant value ρ ∞ corresponding to one of the roots of the equation M(ρ(t))ρ = 0 (42b) ( it is assumed that system (42) has no other kinds of attractors ). Nevertheless, it is instructive to find the asymptotics of the function η(t, τ ) for large t.
For the sake of simplicity, let us take ρ 0 = ρ ∞ , hence ρ(t) = ρ ∞ for any t. The solution of the above problem is sought in the form Then the contribution of the reaction term vanishes, and equation (41b) is reduced to Noting that and using (35a) in Laplace transform of (44) with respect to t ( tilde denotes transformed quantities ), Inverting the asymptotics in (46a), an expression for η is obtained: Note that the distribution function η(t, τ ) is t-dependent for any t, i .e., no stationary distribution over the age τ is established at large t, even when the components densities are constant in space and time. As the time elapses, the weight of the "old" particles grows. The absence of a stationary distribution at large t has important consequences. In particular, it is shown below that the front between two locally stable phases never acquires a stationary shape or a constant velocity. The same conclusion was established formerly for fronts between a stable and an unstable phase [18].

3.3.
Spatially non-uniform problem. For a non-uniform system the general solution of (39a) is Substituting into the definition of ρ (39b) and differentiating with respect to t, we find the following equation: A Fourier transform of (39c) with respect to x yields ( transformed quantities denoted by a hat ) Using the form of the matrix m with σ being the second moment of the step length probability function and D -a diffusion coefficients matrix, the initial condition in Fourier space takes the form or after inversion at the diffusion limit Thus equation (48) becomes Substituting the solution for η(x, t, τ ) into (48) and taking into account the initial condition, an equation for the density of arriving particles η(x, t, 0) reads (53) Denoting y(x, t) = η(x, t, 0) and z(x, t) = t 0 W (τ )η(x, t, τ )dτ , a closed formulation of the problem is presented, including functions of two variables only: 3.4. Approximate solutions. The investigation of the general properties of system (54) is beyond the scope of the present paper. In the rest of the paper the function η(x, t, τ ) is replaced by a trial function, constructed according to the following arguments. Let us consider a front between two large domains where the vector ρ is almost constant, ρ = ρ ± , with ρ ± being two different roots of the equation (42b). In each of these domains η(x, t, τ ) is determined by the relation (46b).

In the intermediate region that function satisfies the condition
In what follows the expression is used as a reasonable approximation that satisfies the aforementioned conditions. Under that assumption it is possible to compute explicitly the integral term in (52): Then (52) takes the form of a normal reaction -diffusion equation with effective diffusion coefficients decaying algebraically with time: The decrease of the effective diffusion coefficients is caused by the "aging" of the ensemble of molecules.
Hereinafter two examples are treated.
3.4.1. Model with time dependent diffusion coefficient and piecewise linear kinetics. Consider a reaction -diffusion equation with a time fractional power diffusion coefficient and piecewise linear kinetics: Such a problem with no reaction was studied in the past [12]. An analytical solution of the problem can be obtained in the case a = 1/2, when the front between the phases ρ = 0 and ρ = 1 is motionless. Defining u = ρ − 1/2 equation (59) is rewritten as Impose an initial condition u(x, t 0 ) = u 0 (x), t 0 > 0 with u 0 (x) satisfying Then u(x, t) possesses the same properties for t > t 0 , enabling one to replace sgn u by sgn x and find an exact solution of (60a). Upon Fourier transform Solving the ordinary differential equation with a general initial condition u(q, 0) = u 0 (q), Hereinafter a simple initial condition is chosen, sufficient to yield a front solution and analyse the effect of anomaly: Inverting the transform and applying (62b), e kτ +q 2 τ γ /γ dτ dq

YANA NEC, VLADIMIR A VOLPERT AND ALEXANDER A NEPOMNYASHCHY
wherein the following identity was used: Note that the erf form includes a singularity at τ = t and a restriction t 0 < t, though it allows to show the monotonicity of the solution ( excluding the point x = 0 ): The limiting values at infinity are ( to leading order ) The steepness of the front is 3.4.2. Model with time dependent diffusion coefficient and cubic non-linearity. As the second example, let us consider a single component model with a cubic nonlinearity ( with the diffusion coefficient scaled to unity ), where u = 2ρ − 1. An analytical solution cannot be presented in that case, but some self-similar solutions can be found in asymptotic limits. An odd solution is considered with the front centre ( the point where u = 0 ) located at x = 0. Let with the parameter α of the similarity variable to be determined. Then First, let us consider the vicinity of the front x = O(1) in the limit t ≫ 1. In that case one finds a self-similar solution of equation (68) i.e. the shape of the front is typical for the Allen-Cahn equation, but its characteristic width decreases with time. Let us consider now the time evolution of the solution in the far tail x ≫ 1 where g = h − 1 is small and (70) can be rewritten as A self-similar solution can be obtained with α = −γ/2 if the reaction term is negligible, hence ( the region where that assumption is justified, is given below ). Then with c being a constant. The estimate of the reaction term (75a) compared with the diffusion or time evolution term shows that the self-similar solution is justified in the region ζ ≫ t 1/2 . Therefore the far tail asymptotics is given by If µ = 0, the similarity solution (69) loses validity, because the centre of the domain wall is not located at x = 0, but moves according to a definite law x = x 0 (t). Instead of (69) let us introduce the ansatz where the function x 0 (t) should be found from the existence condition of the selfsimilar solution. Substituting (78) into (77), an equation similar to (70) is obtained: Similarly to the case µ = 0, a self-similar solution is possible only in the limit t ≫ 1, when the first term on the left-hand side of (79) is negligible. Another ensuing demand The value of v at small µ can be found by means of an asymptotic expansion The leading order equation is giving At the next order an inhomogeneous linear problem is obtained: Thus, the domain wall motion slows down, matching the conclusion obtained for a front between stable and unstable phases [18].

Conclusion.
Existence of front-like solutions was shown for three reactionsub-diffusion models. For the fractional model with piecewise linear kinetics an exact solution for the front shape and an explicit dependence of the front speed on kinetic parameters were obtained. The main observed effect of sub-diffusion is a threshold phenomenon for an anomaly exponent γ close to zero. Specifically, there exists a critical value of the kinetics parameter a = a cr that separates extremely fast ( a > a cr ) and extremely slow ( a < a cr ) fronts. For the fractional model with cubic non-linearity and a small external field the front tails were shown to be asymmetrical due to the memory. The part of the line where the front has not yet passed exhibits a normal exponential behaviour, whereas behind the front the decay is algebraical. The front velocity also changes with γ.
For the integral model based on the continuous time random walk with aging the front was obtained as a similarity solution with a different anomaly dependent scaling for the inner part and the tail. In the presence of an external field the front velocity was shown to decay algebraically with time.
The operators in Fourier space are computed as follows. For a Caputo derivative, defined for the class of differentiable functions, For a Riemann-Liouville type operator and Weyl derivative there is no differentiability condition: Note that the identical operators in Fourier space are obtained under somewhat different demands with respect to the decay at infinity of the function g(t): for Caputo and Riemann-Liouville type derivatives Fourier transform of g(t) should exist, whereas in the case of Weyl derivative the limit in (A.2c ) should also vanish.
Appendix B. Non-existence of complex roots. To show that no complex roots of (11) exist substitute Q ⋆ = re iϕ and separate into two real equations |c| γ r γ cos(γϕ) + k − D r 2 cos(2ϕ) = 0 (B.1a) If c > 0, i.e. the contour is in the upper half plane, for the multi-valued function (−icq) γ it is convenient to choose the branch cut on the negative part of the imaginary axis. Then out of the range [− π 2 , 3π 2 ) the relevant poles might lie within 0 arg q| c>0 π, and since Q = i sgn c q, the argument of Q ⋆ should satisfy − π 2 ϕ π 2 . If c > 0, the branch cut is along the positive part of the imaginary axis, and out of the range [− 3π 2 , π 2 ) the relevant sector is −π arg q| c<0 0 and again − π 2 ϕ π 2 . Therefore for both signs of c the sought roots satisfy |ϕ| π 2 . For a complex root sin(γϕ) = 0 for any 0 < γ < 1, wherefrom it follows also that sin(2ϕ) = 0. Then dividing by sin(γϕ), substituting into (B.1a) and rearranging, D r 2 sin ((2 − γ)ϕ) + k sin(γϕ) = 0. (B.2) Since for any ϕ in the range at hand cos(γϕ) > 0, equation (B.1a) has a solution only if cos(2ϕ) > 0. Then 0 |ϕ| < π 4 , and both terms in (B.2) are either positive, negative or zero, proving that no complex roots Q ⋆ exist.
Appendix C. Algebraical decay of tails. For sgn c = sgn ξ the integration variable in (13) is changed to y = |ξ|q and the integral is split as follows: and all remaining series terms are virtually computed before the limit |ξ| → ∞ is taken, so that the integrals are convergent. In the complementing integral the integration variable is changed y → 1/t: (C.5)