Dynamical features of reaction-diffusion fronts in fractals

Departamento de Medicina, Universitat Internacional de Catalunya, c./Gomera s/n, 08190-Sant Cugat del Valle ́s, Barc lona, Spain Grup de Fı ́sica Estadı ́stica, Departamento de Fı ́sica, Universitat Auto `noma de Barcelona, E-08193 Bellaterrra, Spain Departamento de Fı ́sica, Universitat de Girona, Campus de Montilivi, 17071 Girona, Catalonia, Spain ~Received 12 September 2003; published 29 January 2004 !


I. INTRODUCTION
Reaction-diffusion models are based in general on the partial differential equation ‫ץ‬ P͑x,t ͒ ‫ץ‬t ϭD ‫ץ‬ 2 P͑x,t ͒ ‫ץ‬x 2 ϩ f ͑ P ͒, ͑1͒ where P(x,t) is the density of particles at time t at the position x, D is the diffusion coefficient, and f ( P) is the growth ͑reaction͒ function.The characteristic wave front solutions that arise from these models make them suitable for many different applications which cover biological ͓1͔ and human ͓2͔ invasions, forest fires ͓3͔, epidemics ͓4͔, tumor growth ͓5͔, etc.In spite of this, reaction-diffusion equations are often criticized, by arguing that such simple models cannot account for the complexity of real systems.Specifically, how the intricate features of spatial systems must be modeled is still a current problem.Some efforts have been made in the last years both to show the need for models able to predict spatial complexity ͓6͔ and propose some possible solutions ͓7͔.Among these proposals, one of the most attractive and recurrent ones consists of assuming that the spatial structures exhibit self-similarity properties at a certain range of scales, so fractal scaling formalism may be considered.
In this work, we try to show some of the main concepts involved in adapting reaction-diffusion to self-similar spatial systems.It is well known that the mean-square displacement of a random walker in a fractal object fulfills the subdiffusive behavior ͓8͔ ͗x 2 ͘ϳt 2/d w ,

͑2͒
where d w у2 is the random-walk dimension of the fractal ͑for d w ϭ2 one recovers the classical case͒.The main advances in the field of transport in fractals have been focused on finding a transport equation for the probability density P(x,t) of finding a random walker at time t at a distance x from its starting point.At present, it is known that for all particles starting from the same origin and for sufficiently large times and distances, the following form of the probability density ͓9,10͔ is valid for a large class of fractals.In Eq. ͑4͒, d min is the fractal dimension of the minimum distance between points of the fractal ͓11͔.The classical, Euclidean case corresponds to d w ϭ2,d min ϭ1 and, therefore, uϭ2 so that the solution is a Gaussian.
Transport equations proposed to date have tried to reproduce results ͑2͒-͑4͒ somehow.The first attempt was made by O'Shaugnessy and Procaccia ͑SP͒ ͓12͔.They derived from scaling and renormalization arguments the SP equation

͑5͒
where d f is the geometric fractal dimension of the fractal, D(x)ϭD*x 2Ϫd w , and the constant D* is a kind of diffusion coefficient.The exact solution for this equation is known and its second moment behaves like t 2/d w , in agreement with Eq. ͑2͒.However, this solution lacks the scaling ͑3͒ and ͑4͒ predicted by numerical simulations.Later on, Giona and Roman ͑GR͒ constructed the fractional GR equation ͓13͔ To overcome these difficulties, more recently a generalization of the SP equation including a temporal fractional derivative Finally, very recently we have proposed a partial differential equation ͓the Campos-Me ´ndez-Fort ͑CMF͒ equation͔ from a stronger physical justification than the previous models, which also improves the results obtained by previous approaches ͓9͔.
Our goal is to find the analytical relationship for the speed of fronts when the above transport equations couple to a reaction process modeled by logistic-KPP ͑Kolmogorov-Petrovskii-Piskunov͒ kinetics.We employ the method of reduction of the reaction-transport equation to a Hamilton-Jacobi.Also, we are able to derive in some cases the speed of the fronts by using the local-equilibrium ͑LE͒ approach ͓15͔ when the spatial correlations are small.

II. SP EQUATION WITH REACTION
First of all we consider that the reaction-diffusion process in a fractal is described by the equation for the probability density of O'Shaughnessy and Procaccia ͓Eq.͑5͔͒ coupled to a KPP kinetic term ‫ץ‬ ‫ץ‬x ͩ D*x d f Ϫd w ϩ1 ‫ץ‬ P ‫ץ‬x ͪ ϩrP͑1Ϫ P ͒. ͑8͒

A. Reduction to a Hamilton-Jacobi equation
In order to find the asymptotic speed for the traveling wave fronts in Eq. ͑8͒ we will first make use of Hamilton-Jacobi dynamics ͓16͔.The starting point is the hyperbolic scaling procedure t→t/, x→x/, and the representation of the rescaled probability density function P (x,t) ϭ P(x/,t/) in WKB form where the action functional G has to be found.
The first and third terms in the right-hand side of Eq. ͑10͒ have the same order of magnitude and in the asymptotic limit (→0) both terms may be neglected in front of the second term, and in consequence the Hamilton-Jacobi for the front propagation in a fractal media is where G(x,t)ϭlim →0 G (x,t).In order to find the solution for G(x,t) in Eq. ͑11͒ we make use of the Hamilton equations where HϵϪ‫ץ‬G/‫ץ‬t and p()ϵ‫ץ‬G/‫ץ‬x() and stands for the temporal coordinate.The solution of the Hamilton-Jacobi equation ͑29͒ can be written as where L(x,)ϭp()͓dx()/d͔ϪH is the Lagrangian associated with H. Integrating Eq. ͑12͒, one has x() ϭx(/t) 2/d w under the boundary conditions x(0)ϭ0, x(t) ϭx, and and from Eq. ͑13͒ The speed of the front is computed from G(x,t)ϭ0, and after inverting the hyperbolic scaling one finally has the exact expression for the speed of the front which describes a decelerated front, and we note that the acceleration depends on the parameter d w .

B. Local-equilibrium approach
The local-equilibrium hypothesis has been very successful in thermodynamics ͑where one assumes a Gibbs equation with the local values of temperature, pressure, etc.͒, radiative transfer, astrophysics ͑where one assumes a radiative emission-absorption equilibrium at the local temperature͒, etc.Here, we will apply such a heuristic approximation to the problem of front propagation.Equation ͑8͒ can be rewritten as This equation has the form which is an advection-reaction-diffusion equation in a medium moving at speed V(x).For homogeneous values of D and V, the speed of fronts is obviously vϭ2 ͱrDϩV,

͑17͒
which is nothing but Fisher's speed 2ͱrD ͑which holds for an observer attached to the moving medium͒, as seen by an observer moving with speed ϪV.This result vϭ2ͱrDϩV may be also derived by using the linearization and variational techniques in Ref. ͓17͔.We propose that, if inhomogeneities are sufficiently smooth, this equation should hold locally, so that we obtain the local-equilibrium prediction which for Eq.͑15͒ yields Note that if one takes the hyperbolic scaling x→x/, t →t/ with →0 in Eq. ͑19͒, one observes that the second term on the right-hand side ͑rhs͒ of ͑19͒ is negligible, compared to the first one, and Eq.͑19͒ can be written as ͱrD* x (d w Ϫ2)/2 for x,t→ϱ, which may be integrated to yield x(t)ϭ͓d w ͱrD*t͔ 2/d w , where we have assumed the initial condition x(0)ϭ0.Finally, the speed may be obtained as

͑20͒
which coincides with Eq. ͑14͒.The asymptotic speed does not depend on the fractal dimension d f , which could be something surprising at first sight.In order to further explore this point, we have solved numerically Eq. ͑8͒, under the initial condition P(x,0)ϭ1 for xϽ0 and P(x,0)ϭ0 for x Ͼ0, for two very different values of d f and we have observed, as we show in Fig. 1, that the speed depends weakly on d f only at the very initial transient and loses this dependence as time grows.Moreover, in Fig. 1 we compare the speed given in Eq. ͑14͒ or ͑20͒ with the speed obtained from Eq. ͑19͒.
To do this we must first solve Eq. ͑19͒.Under x(0)ϭ0 we obtain tϭ y where xϭy 2/d w .In Fig. 1 we have solved numerically the transcendent equation ͑21͒ and evaluated dx/dt to obtain the speed.Note that the effect of d f ͓which appears in Eq. ͑21͒ but not in Eq. ͑20͔͒ on the asymptotic speed of the front is inappreciable.This shows that it is justified to neglect the second term in the rhs of Eq. ͑19͒, as done above.

III. GR EQUATION WITH REACTION
In this section we study the speed of fronts for the fractional advection equation GR ͑6͒ with reaction where FIG. 1.Comparison between numerical results for the speed of fronts of Eq. ͑8͒ for d f ϭ1 and for d f ϭ0.01 and the analytical results from the HJ method given in Eq. ͑14͒ ͑solid line͒ and the numerical result for the LE method obtained from Eq. ͑21͒ with d f ϭ1.All magnitudes are dimensionless.We observe a very good agrement between analytical and numerical solutions in the asymptotic regime.It is clear that d f does not affect the value of the speed for large times.We have taken D*ϭrϭ1 and d w ϭ2.32.
and zϭtϪtЈ.As in the preceding section we first take the asymptotic limit for large space and time by employing the hyperbolic scaling.The left-hand side of Eq. ͑6͒ is where P (x,0)ϭ0 for xϾ0, and making use of Eq. ͑9͒ one finds The right-hand side of Eq. ͑6͒ is transformed into Taking →0 one gets the Hamilton-Jacobi equation From the Hamilton equations one has the Lagrangian and the speed of the front will be , which is time independent.

IV. FRACTIONAL SP EQUATION WITH REACTION
Another equation proposed is the fractional version of the SP equation ͑7͒ which with reaction takes the form where the fractional time derivative is defined as ͓14͔ Taking the hyperbolic scaling and Eq.͑9͒ in the limit →0 one finds the Hamilton-Jacobi equation From the Hamilton equations one obtains p() ϭBx() Ϫ1ϩd w /2 and where B is an integration constant to be determined from the condition x(ϭt)ϭx, and therefore On the other hand, the Lagrangian function is From G(x,t)ϭ0 we have and inserting this into Eq.͑25͒ one has and taking dx/dt the speed is

͑26͒
once the hyperbolic scaling is inverted.Note that for ␥ϭ1, Eq. ͑14͒ is recovered.However, the scaling law for the speed of the front in time is the same as for ␥ϭ1 in Eq. ͑20͒ so that the fractional derivative, although affecting the mean-square displacement, does not affect this scaling law.

V. CMF EQUATION WITH REACTION
We start from the CMF equation

͑27͒
which has been derived very recently ͓9͔.

A. Reduction to a Hamilton-Jacobi equation
After taking into account the hyperbolic scaling and the field G (x,t)ϭϪ ln P(x/,t/) one has from Eq. ͑27͒

‫ץ‬G
‫ץ‬t The first term in the left-hand side ͑lhs͒ of Eq. ͑28͒ and the second term in the rhs of Eq. ͑28͒ are O( 0 ).The second term in the lhs and the first term in the rhs of Eq. ͑28͒ are O( ␣ϩ1 ) while the third term in the lhs is O( ␣ ).For fractals d min у1, so that uуd w /(d w Ϫ1) or ␣Ͼ0.Therefore, in the limit →0 one has G(x,t)ϭlim →0 G (x,t), lim →0 f (e ϪG / )ϭ0, provided G (x,t)Ͼ0 and By keeping the terms up to O( ␣ ), one has

͑29͒
Equation ͑29͒ is the Hamilton-Jacobi equation for the problem ͑27͒ and may be solved by using Hamilton's equations Integrating Eq. ͑30͒, one has x()ϭx(/t) 2/d w under the boundary conditions x(0)ϭ0,x(t)ϭx, and Finally, from Eq. ͑13͒ we obtain If one inverts the hyperbolic scaling in Eq. ͑31͒, by taking x→x and t→t, one has The speed of the front is computed from G(x,t)ϭ0, and after inverting the hyperbolic scaling one finally has It is important to note that the scaling law in Eq. ͑32͒ does not depend on d f and d w but only on d min

B. Local-equilibrium approach
Following the same method as in previous sections, we find for Eq.͑27͒

͑33͒
which leads us to the front speed ϳt (1/d min )Ϫ1 for t→ϱ.͑34͒ It is interesting to note that although the LE prediction ͑34͒ has the same scaling law as the HJ result ͑32͒, the factor is different unless d w ϭ2d min .However, we have checked numerically that for the typical ranges of d w ͑from 2 to 4) and d min ͑from 1 to 2) in fractals, the differences between the two models are negligible ͑see Fig. 2͒.

VI. DISCUSSION
After all the mathematical formalism, the main conclusion we obtain is that, despite all these equations seeking to describe the features of transport on fractals, the characteristics of the fronts predicted in the four cases shown are clearly different.First of all, we should note that a diffusion equation on fractals must agree with the results ͑2͒-͑4͒.We have shown ͓9͔ previously that equation CMF is the only one that can exactly reproduce them, so we expect that the fronts predicted by this equation correspond to the real case.
One may wonder why fronts on fractals should be accelerated.To answer this we need to introduce in our discussion the ''chemical distance'' l ͓11͔, which is defined as the shortest path between two points belonging to the fractal ͑Fig.3͒.
Usually, random walks ͑or diffusion͒ on a fractal is described as follows.A particle at point s ͑see Fig. 3͒ can jump to every one of its first neighbors with the same probability.At the next time step, it may again jump to the new first neighbors with the same probability, and so on.We can see that the possible points reached after two jumps forward have in common that their chemical distance to the origin s is the same ͑but not their Euclidean distance͒.Hence, we see that transport on fractals take place through the chemical distance space and in this space the fractal fronts are expected to show constant speed, as in homogeneous media.
Nevertheless, we are usually interested in the results for the Euclidean space, so we need the well-known relationship between r and l ͓11͔, where the condition d min у1 comes directly from the fact that l is always greater than r, as seen in Fig. 3. From this expression we can conclude that, assuming that fronts advance at a constant speed in the chemical distance space, they cannot do so in the Euclidean space.Moreover, we observe that this behavior is due to the parameter d min ͑in the nonfractal case we have d min ϭ1, so the behavior in both spaces will then be the same͒; it agrees with the CMF equation ͑27͒, which predicts an acceleration dependence only on d min ͓see Eq. ͑32͒ or ͑34͔͒, while the other equations analyzed ͑Secs.II, III, and IV͒ do not take into account this essential parameter.
In fact, the form of the time exponent in Eqs.͑32͒ and ͑34͒ can be justified.We may define l f r and r f r as the chemical and Euclidean distances of the wave front position, respectively.As the front speed is constant in the chemical distance space, we can assume that l f r grows linearly with t.This relation, in addition to Eq. ͑35͒, leads to r f r ϳt 1/d min , ͑36͒ and the time dependence expected for the fronts in the Euclidean space is then vϳ r f r t ϳt (1/d min )Ϫ1 , ͑37͒ in agreement with the acceleration predicted by the CMF equation ͓Eq.͑32͒ or ͑34͔͒.This result, which had been already predicted for the speed of propagation of fronts in percolation clusters ͓18͔, should be valid for all those transport processes on fractals which, as it happens in most cases, take place through the chemical distance space ͑as in Fig. 3͒.
In conclusion, we have shown that the CMF equation not only describes the features of diffusion on fractals better than previous ones ͓9͔, but it also predicts some essential features of the propagative processes on those heterogeneous media, i.e., the speed and acceleration of the wave fronts derived when a reaction ͑logistic͒ process, widely used in biophysics ͓19͔, is considered.In consequence, we think that CMF equation is the best analytical approach proposed to date for description of transport on fractals.FIG. 2. Comparison between numerical results for the speed of fronts of Eq. ͑27͒ and the analytical results from the HJ method given in Eq. ͑32͒ ͑solid line͒ and the numerical result for the LE method obtained from Eq. ͑34͒ with d f ϭ1.All magnitudes are dimensionless.There is good agrement between analytical and numerical solutions in the asymptotic regime.We have taken D*ϭr ϭ1 and d w ϭ2.32, d min ϭ1.1.
One of the main results presented here is the wave front speed ͑32͒, which shows that the acceleration is determined by the parameter d min , as argued theoretically above.It is interesting to note that this parameter is not the same as that responsible for anomalous diffusion, namely, d w ͓see Eq. ͑2͔͒, contrary to what one could expect.We think that an exhaustive analysis of the meaning of the parameter d w , which has been traditionally defined just from Eq. ͑2͒, is still needed.This and many other questions which have not been explained theoretically yet ͓18͔ show very clearly that there is still a lot of work to do on fractal dynamics.Although theoretical research on fractals decreased after a boom in the early 1980s, we consider that efforts in this field, as the one presented here, are still strongly useful.

FIG. 3 .
FIG.3.Sierpinski gasket as an example of random walk on a fractal.The bold lines show the difference between the Euclidean distance r and the ''chemical distance'' l between two points of the fractal s and sЈ ͑note that lуr for any couple of points͒.