Reaction-diffusion wave fronts on comblike structures

From the Hamilton-Jacobi formalism, an explicit expression for the speed of wave front propagation along the backbone of comblike structures is obtained. This expression, through the waiting-time distribution function, takes into account the number of sites and their distribution in the secondary branches. Our theoretical results are supported by numerical simulations of the reaction random-walk process on the structure. Finally, a more complex situation such as the Peano basin structure is also considered, both theoretically and numerically, exhibiting a good agreement too.


I. INTRODUCTION
Lattice models and some other discrete models have become very popular in recent years for dealing with a wide variety of dynamical systems ͓1,2͔.In this work we want to show an analytical method for studying the propagation of reaction random-walk wave fronts through heterogeneous lattices when these structures are made up by a main backbone with a regular distribution of secondary branches.Thus, our work is especially appropriate and easy to understand for comblike models ͑see Fig. 1͒ and so we focus on this kind of structure, which has become of special interest for many researchers in recent years ͓3͔.
For our purpose, we start from the continuous-time random walk ͑CTRW͒ framework ͓4͔.Although this approach is usually applied to continuous systems, we show that the discrete nature of lattice models can be introduced there by choosing the appropriate distribution functions describing the jumps made by the walkers.By doing that, our models will be able to capture the properties of the underlying random-walk process in the lattice ͑and, in the limit of large distances, we can interpret it as a diffusion process͒.
After that, we introduce in the CTRW a reaction process taking into account the creation ͑annihilation͒ of walkers.This yields, as usual, traveling wave front solutions spreading through the media.These models are analyzed by means of the Hamilton-Jacobi techniques, which in recent years have been proved quite useful for the determination of wave front speeds on reaction random-walk systems ͓5,6͔.The results for the front speeds are compared with random-walk numerical simulations performed within the comb structures, showing in general a good agreement.Finally, we justify that our method can be applied to a large range of lattice models, provided that we are not interested in the propagation through the whole structure but just in the direction of the backbone.

II. RANDOM WALKS
A diffusion process can be expressed in the CTRW framework in terms of the density of particles ͑x , t͒ as where ⌿͑x , t͒ is the probability for a walker to perform a jump of length x when a time t has elapsed since its previous jump.Usually, the distribution of jump lengths and the distribution of time probabilities are independent, so the function ⌿͑x , t͒ can be written in the decoupled form ⌿͑x , t͒ = ͑t͒⌽͑x͒ and Eq.͑1͒ turns into where ͑t͒ represents the distribution function for waitingtime probabilities and ⌽͑x͒ the distribution of jump lengths.We first consider the case of a discrete one-dimensional chain where the first neighbors are separated by a distance ⌬x.Then, a single random walk, where each walker moves to one of its first neighbors with equal probability after a remaining time , is characterized by the distributions

͑4͒
so the time and the space in the model are taken as discretized.In this way, discrete systems and lattices can be analyzed by the CTRW approach.Next, we add to every site of the backbone a secondary branch of length l, so a comblike structure appears ͑Fig.1͒.In this situation, a walker that is at a site on the backbone can spend a certain time in the secondary branch before passing to one of the first neighbors in the backbone.Therefore, if we are only interested in the behavior of the system in the direction of the backbone, then we can interpret the secondary branches as introducing a delay time between the sites in the backbone.In this case, the random walk within the comb structure is given by Eq. ͑4͒ and a new distribution ͑t͒ which has to include the effect of the delay by the secondary branches.
To determine analytically the effect of the branches we can apply some convolution rules which were partially discussed by Van den Broeck ͓8͔ for the case of homogeneous lattices.
͑i͒ Let us consider that a walker is initially at a certain site within the secondary branch.If the walker moves further on the secondary branch ͑so it goes away from the backbone͒ the probability for it to return to the initial site in a time t is a convolution of factors ͑then it becomes a product in the Laplace space͒.
͑ii͒ The total probability for that walker to return to the initial site is determined by the sum for t from 0 to ϱ.
͑iii͒ When the walker reaches a crossing ͑so it can choose between different ways͒ the total probability is the sum of the probabilities for each possible way.
Note that for the case of comb structures like that in Fig. 1, the secondary branches have no crossings, so rule ͑iii͒ is not necessary there; however, we will show in Sec.IV that in many other situations it must be taken into account too.
For the sake of generality, let us consider that when the walker is at a site in the backbone it can jump to another site in the backbone with probability ␣ or go into the secondary branch with probability 1 − ␣.This generalization allows us to analyze, for example, structures like those shown in Sec.IV below as equivalent to that in Fig. 1 just by choosing the suitable value for ␣.
Now we can assume that initially a walker is located in the backbone and apply the rules ͑i͒-͑iii͒ to determine ͑t͒.We will examine three specific cases.
(a) Comb structure with l = ⌬x.In this case there is only one site within the secondary branch, so the walker can only jump in the direction of the backbone with probability ␣ or go into the branch with probability 1 − ␣ and then come back to the initial site at the next jump.In consequence, the time it takes to reach one of the first neighbors in the backbone is t = with probability ␣, t =3 with probability ͑1−␣͒ ϫ 1 ϫ ␣, t =5 with probability ͑1−␣͒ 2 ϫ 1 2 ϫ ␣, and so on.Hence, we can write intuitively the general form ͑t͒ as

͑5͒
The rules shown above for ͑t͒ should now reproduce this behavior.For this purpose, we need to work in the Laplace space and so we will use ˆ͑s͒, which is defined as the Laplace transform of ͑t͒.
The rules ͑i͒-͑iii͒ allow us to write the expression where ˆ0 is the probability distribution for one single jump ͓in our case, ˆ0 = exp͑−s͒, which is the Laplace transform of Eq. ͑3͔͒.
Equation ͑6͒ is derived as follows.The term ͑1−␣͒ ˆ0 2 within the sum represents, according to rule ͑i͒, the probability function for each time the walker goes into the secondary branch.This expression must sum up to infinity ͓rule ͑ii͔͒ to take into account that the walker can go into the branch 1,2, ... ,ϱ times.Finally, the factor ␣ ˆ0 for the final jump to the first neighbor in the backbone is added.
The next step is to calculate the antitransform of ͑6͒ in order to compare it with the predicted result in ͑5͒.Unfortunately, that antitransform is not possible to find analytically.Instead of that, it is easy to see that the expression ͑6͒ may be written ͑by Taylor series͒ as Now it is straightforward to see that the Laplace transform of ͑5͒ is ͑7͒ and so our method for determining ˆ͑s͒ is proved to be valid in this case.
(b) Comb structure with l =2⌬x.Analogously to the case seen above, if now the secondary branch is two sites long we can write the distribution for the time probabilities as

͑8͒
In this equation, a new sum for the index k appears because now the walker can go twice away from the backbone and for each one of them we must apply rule ͑i͒.We have also used the fact that within the linear secondary branch the jumps to the first neighbor have probability 1 / 2.
We have checked the result in Eq. ͑8͒ by random-walk simulations within a comb structure.According to the results so obtained ͑not shown͒, the antitransform of expression ͑8͒ agrees exactly with the jump probabilities to the first neighbor as a function of time, and so it confirms the validity of our method.
(c) Comb structure with l → ϱ.We have mentioned that, as long as the walker moves away from the backbone, a new convolution factor appears in ˆ͑s͒.So that, for the case l → ϱ, we would have theoretically infinite convolution factors in the expression for ˆ͑s͒.
Nevertheless, we can greatly simplify this situation.Consider that the walker is at the first site in the secondary branch ͑point AЈ in Fig. 1͒ and moves away from the backbone; then, we can call A Ј the probability distribution of returning for the first time to the point AЈ after a time t.Now imagine the same situation but for the initial point AЉ; it is easy to see that, as l → ϱ, the condition A Љ → A Ј has to hold.
Thus it is possible to use the rules ͑i͒-͑iii͒ to determine A Ј .By doing this, we obtain the expression which is an expression equivalent to the form ͑6͒ with ␣ =1/2 ͑because within the secondary branch every jump to a first neighbor has probability 1 / 2͒.Now we can introduce the condition A Љ = A Ј , which is strictly correct for l = ϱ.
Solving so Eq.͑9͒ yields Once we have found this result, the distribution ˆ͑s͒ comes straightforwardly by applying again the three rules ͑i͒-͑iii͒ above.It leads us to

͑11͒
an expression which has been checked by simulations too, confirming that Eq. ͑11͒ agrees exactly with the results expected.
At the sight of this agreement found in the cases ͑a͒-͑c͒, we can affirm that our method is able to describe exactly the random-walk process in these structures along the backbone.

III. REACTION-RANDOM WALKS
Our main objective, as said above, is to use the probability distributions we have found to determine the properties of reaction random-walk processes and the wave front solutions they yield.In order to extend the CTRW approach shown above to the case of reaction random-walk systems we will add to Eq. ͑2͒ a term that accounts for the reaction process.So we obtain the mesoscopic equation where the function ͑t͒ is defined as the probability to remain still at least a time t before performing a new jump, In consequence, the reaction function f, which determines the characteristics of the reaction process, is applied for all the particles that are at position x and arrived there a time tЈ ago.For a more rigorous description, the reader can find in ͓6͔ how Eq. ͑12͒ can be derived straight from the master equations in the CTRW framework.
To analyze the reaction random-walk equation we first need to give a specific form to f.A reasonable choice which is not too difficult to deal with and has been used before successfully in many different fields, such as population dynamics, tumor growth ͓9͔, and so on, is the Fisher-Kolmogorov-Petrovskii-Piskunov equation where a is a constant growth rate.
The system ͑12͒-͑14͒ is expected to yield in general wave front solutions which connect the unstable ͑ =0͒ and the stable ͑ =1͒ state in Eq. ͑14͒.Thus, the speed of the wave front becomes the essential parameter to find in order to determine the behavior of the system.
Recently, we have shown ͓6,10͔ that the wave front speed for the case ͑12͒-͑14͒ can be obtained analytically by Hamilton-Jacobi techniques ͓5͔; as a result, we can write the expression for that speed as where H͑p͒ or p͑H͒ is the solution of the Hamilton-Jacobi equation

͑16͒
The distributions ˆ͑H͒ and ⌽ ˆ͑p͒ are defined as so they are the Laplace transform and the bilateral transform of the waiting time and the jump length distributions, respectively.
The importance of the result ͑15͒-͑17͒ lies especially in its generality.We can introduce any pair of distributions ͑t͒ , ⌽͑x͒ determining a diffusion pattern and obtain as a result the speed for the corresponding reaction-diffusion fronts.Now, we want to apply the Hamilton-Jacobi result for random walks in the comb structures studied in Sec.II.According to our initial hypothesis we only consider the random-walk process through the backbone and interpret the secondary branches as included into the definition of ͑t͒.Thus we can take the jump length distribution ͑4͒ and the time distributions obtained in Sec.II ͓note that, as we worked there in the Laplace space, we already know the form of ˆ͑H͔͒.
Introducing Eq. ͑4͒ into Eq.͑16͒ and solving for p we obtain from the first equation in ͑15͒ the front speed for the three cases in Sec.II, with ˆ͑H͒ = ˆ͑s = H͒ given by Eqs.͑6͒, ͑8͒, and ͑11͒.In Eq. ͑18͒, the minimum that gives us the value for v cannot be computed analytically; instead of that, we must find it by numerical methods.This operation has been performed in Fig. 2, where we compare the results found with the speed obtained from simulations within comb structures where we have implemented a reaction random-walk process.The details for the simulations are given in ͓7͔; basically, we simulate a random-walk process within the structure and then the reaction function ͑14͒ is applied to every site of the structure at every time step of the simulation in order to implement the reaction process.From Fig. 2, we can conclude that the agreement between Eq. ͑18͒ and simulations is excellent for low values of a.Nevertheless, as a grows there appear some little discrepancies between the two results ͑for high values of a the simulations become more difficult and it prevents us from extending the results in this regime͒.
These differences may be explained as follows.For the random-walk process, we have proved in Sec.II that our method agrees exactly with simulations.On the contrary, the reaction process is implemented differently in both cases: in the simulations we apply the reaction function to every site in the lattice, but our theoretical method considers that all the walkers are concentrated in the backbone.So, when we apply f͑͒, the dependence in introduce differences between the two situations.However, we see in Fig. 2 that only if the reaction process is important ͑for high values of a͒ do the differences between the two situations become apparent.This deviation increases with l and it is maximum for l → ϱ because for this case there are more sites in the structure that do not belong to the backbone.See Fig. 3 Finally, we note that in the regime for high values of a the wave front speed tends asymptotically to the speed of the individual particles ⌬x / ͑which is equivalent to v / ⌬x → 1 in Fig. 2͒, a behavior that we discussed in ͓6͔ and that now we have proved analytically ͑see the Appendix͒.

IV. OTHER CASES
Here, we have focused on the study of wave front solutions across comblike structures because they are specially suitable to illustrate our method.However, we want to point out that we are allowed to deal with a wide variety of heterogeneous lattices.The one important restriction is that we must be only interested in the dynamic behavior across a backbone and then consider the rest of the structure as secondary.
One of the most interesting cases we could analyze ͑especially for practical purposes͒ would be fractal structures.In Fig. 3 we show the construction for the well-known Peano basin ͑as a function of the construction level Q͒, which has been used before for the modelization of fractal river basins ͓11,12͔.For instance, if we are only interested in the behavior in the direction of the backbone, then the case Q = 2 for the Peano basin leads us to the waiting-time distribution function

͑19͒
where we have used ␣ =1/2 ͑the same jump probability in all directions͒.The wave front speed derived from Eq. ͑19͒ can be obtained again by means of Eq. ͑18͒.In the plot in Fig. 4, we show these speeds for the case Q = 2 and also for the case Q =5 ͓ ˆ͑s͒ is not shown here for simplicity͔, the highest level we have analyzed by our method.We compare them with reaction random-walk simulations within a Peano basin with Q =10 ͑for simulations, we need a high value of Q to allow the formation of fronts; although the value of Q is different in the theory and the simulations, we observe that the results for v converge very fast and so taking the greatest values for Q in the theory would not improve our results substantially͒.In consequence, we find for the Peano basin the same behavior as that observed for comb structures ͑that is, a good agreement between theory and simulations for low values of a and some discrepancies as a grows͒.
We also mention that percolation clusters, one of the bestknown fractal structures, are another case for which it has been proposed that comblike models can be a good approach to study their properties ͓13͔, though in this case there exists some controversy about it.

V. CONCLUSIONS
From the CTRW framework we have used a mesoscopic equation to describe the front propagation in comblike structures in the backbone direction.We have exposed how to derive the waiting-time distribution function in the Laplace space, which has to take into account the number of sites in the secondary branches.These are assumed to emerge from the backbone's sites separated by a distance ⌬x.The Peano basin and three particular cases of comblike structures have been studied.We have employed the Hamilton-Jacobi formalism to obtain an explicit expression for the speed of front propagation in these structures.Our theoretical results have been compared to numerical simulations exhibiting a good agreement.The structures we have analyzed are just some examples to show the wide range of application for our method.We think there are still many more different structures and lattices with practical interest and whose dynamic properties can be successfully analyzed by the approach we have presented here.In Ref. ͓6͔, we argued that there is a physical restriction to the speed of reaction-random-walk fronts.This restriction was called the finite jump speed effect and it says that the speed of the front cannot be higher than the speed of the individual particles in the process.The speed of the front tends to increase as the reaction term gets higher, so in our case, where the reaction term is controlled by the value of a, we expect that in the regime a → ϱ the speed of the fronts will tend asymptotically to the speed of individual particles.Next, we shall prove it analytically.
We observe from the expression ͑18͒ for the speed of fronts that a → ϱ implies H → ϱ too ͑from the condition H Ͼ a͒.So, first of all we must look for the expression of ͑t͒ in this limit of high H.
For this purpose, we can observe that a → ϱ means that the reaction process is almost immediate ͑the characteristic time for the process is a −1 → 0͒, so in this regime the reaction process is dominated by the first particles jumping, that is, those particles that jump after a waiting time t = .According to our arguments in Sec.II, this means that we can write where ␣e −H represents the distribution function corresponding to the probability ␣ of jumping to the first neighbor of the backbone after a time .It can be checked that Eq. ͑A1͒ holds in all the cases reported in this article-Eqs.͑6͒, ͑8͒, ͑11͒, and ͑19͒.From Eq. ͑A1͒ and the condition a Ͻ H we can also write Now it is easy to see that H → ϱ leads us to v → ⌬x / , so the front speed tends in this limit to the value of the individual jump speed.

FIG. 1 .
FIG.1.Representation of a comb structure with distance between first neighbors ⌬x and secondary branches of length l.The symbols A , AЈ , AЉ , . . .are the names we give to some sites of the lattice ͑see text͒.

FIG. 2 .
FIG.2.Comparison of the wave front speeds ͑normalized by v / ⌬x͒ found from the Hamilton-Jacobi method ͑lines͒ and from the simulations within the comb structures ͑circles͒ as a function of the reaction parameter a for the three cases studied here.All the parameters in the plot are dimensionless.

ACKNOWLEDGMENTS
Daniel Campos acknowledges the Departament d'Universitats, Recerca i Societat de la Informació of the Generalitat de Catalunya.This work was partially funded by the Generalitat de Catalunya under Grant No. SGR-2001-00186, and by the MICYT under Grants No. REN-2003-00185 CLI and No. BFM 2003-06033.APPENDIX: WAVE FRONT SPEED IN THE REGIME a \ ؕ using Eqs.͑A2͒ and ͑A3͒, we can write the speed v from Eq.

FIG. 4 .
FIG. 4. Wave front speeds ͑normalized by v / ⌬x͒ for the Peano basin, analogously to Fig. 2. Circles represent values from simulations with Q = 10 and lines are obtained from Eqs. ͑18͒ and ͑19͒ for the cases Q = 2 and 5.All the parameters in the plot are dimensionless.