Fronts from complex two-dimensional dispersal kernels : Theory and application to Reid ’ s paradox

Bimodal dispersal probability distributions with characteristic distances differing by several orders of magnitude have been derived and favorably compared to observations by Nathan et al. Nature London 418, 409 2002 . For such bimodal kernels, we show that two-dimensional molecular dynamics computer simulations are unable to yield accurate front speeds. Analytically, the usual continuous-space random walks CSRWs are applied to two dimensions. We also introduce discrete-space random walks and use them to check the CSRW results because of the inefficiency of the numerical simulations . The physical results reported are shown to predict front speeds high enough to possibly explain Reid’s paradox of rapid tree migration. We also show that, for a time-ordered evolution equation, fronts are always slower in two dimensions than in one dimension and that this difference is important both for unimodal and for bimodal kernels. © 2007 American Institute of Physics. DOI: 10.1063/1.2733631


I. INTRODUCTION
Reaction-diffusion phenomena are observed in many applied biophysics systems, e.g., tumor growth, 1 patterning in microchannel flows, 2 nanobiosensors, 3 etc.However, the diffusion ͑or Laplacian͒ approximation does not hold in some cases-not only in matter transport but also in heat conduction 4 and radiative transfer. 5,6Indeed, classical diffusion ͑Fick's law͒ is well known to break down in some relevant phenomena.One example is time-delayed diffusion due to the rest time of particles ͑or individuals͒ between successive jumps 8,9 ͑this delay can be the same for all jumps [10][11][12] or there may be several possible delays [13][14][15] ͒.A second case is sometimes called reaction dispersal.Broadly speaking, reaction dispersal refers to dispersal probability distributions ͑kernels͒ which do not vanish for very long jump distances ͑compared to the front width͒.This leads to the break down of the diffusion ͑or second-order͒ approximation.7][18][19] An especially important subcase is that of kernels with several components with characteristic distances differing by several orders of magnitude.It has been long suspected that such kernels may explain a very important, unsolved biophysical problem, namely, the fact that the observed speeds of forest postglacial recolonization fronts are much faster than those predicted by single-kernel reaction-dispersal models.This disagreement is called Reid's paradox. 18Many authors have shown that hypothetical longdistance dispersal ͑LDD͒ events could solve Reid's paradox, 18 but all previous approaches are based on kernels fitted to short-distance data, with purely hypothetical LDD events 18,19 ͓and almost always using one-dimensional ͑1D͒ models͔.However, Nathan and co-workers have derived and tested very interesting kernels with two components: a shortdistance component ͑of the order of 10 m͒ and a second, very rare, LDD component ͑covering distances of 10 3 -10 4 m but observed for only about 0.2% of seeds released from the parent tree͒.They derived such bimodal kernels by means of fluid dynamics simulations of atmospheric transport including turbulent-uplifting events that had been previously neglected.They also checked their new kernels by comparing predicted vertical deposition patterns and uplifting probabilities to observed data. 7,20,21This opens the possibility to explain Reid's paradox using complex dispersal kernels which are not hypothetical but derived from physical principles.Below we derive the appropriate front speed formulas for such complex kernels.We also show that the predicted front speeds are about 10 2 -10 3 m/yr ͑which are two orders of magnitude higher than those obtained neglecting the LDD component͒.This could eventually solve Reid's paradox.

II. EVOLUTION EQUATION
Many two-dimensional ͑2D͒ reaction-dispersal models are based on the equation [8][9][10][11][12]22 p͑x,y,t where p͑x , y , t͒ is the population number density at position ͑x , y͒ and time t.The dispersal kernel ͑⌬ x , ⌬ y ͒ is the probability per unit area that a particle ͑or individual͒ that was at ͑x + ⌬ x , y + ⌬ y , t͒ jumps to ͑x , y , t + T͒.T is the time interval between two subsequent jumps ͑usually T =1 generation 10,11 ͒.The last term, R͓p͑x , y , t͔͒, corresponds to the reaction ͑or reproduction͒ process.In contrast to, e.g., human populations, 10,14 dispersal ͑of seeds͒ for tree populations takes place during a specific period of the year only, always after reproduction ͑seed production͒.The dispersal and reproduction are not independent, so a time-ordered evolution equation for the adult tree number density must be a͒ Electronic mail: joaquim.fort@udg.esused, instead of Eq. ͑1͒.The time-ordered evolution equation is 16,17 p͑x,y,t where R 0 is the net reproductive rate ͑number of seeds per parent tree and year which survive into an adult tree͒, and T is the age ͑in yr͒ at first reproduction ͑T is sometimes called the generation time͒.Equation ͑2͒ is exact only for species with nonoverlapping generations ͑i.e., such that parents reproduce only once and then die͒. 16,17But previous 1D results with hypothetical kernels show that substantially more complicated, age-structured models do not change the order of magnitude of the front speed. 18,19Thus, the interesting case of age-structured extensions of the 2D Eq. ͑2͒ will be tackled in future work.Here we will use the approximate Eq. ͑2͒.
This will make it much easier to focus our attention on the main physical question posed in the Introduction: How do the dispersion kernels in Ref. 7 affect the speed of fronts?
The kernels in Ref. 7 have two components: a short-distance one, "S" ͑with probability p S and distances ϳ10 2 m͒, and a long-distance one, "L" ͑covering distances of 10 3 -10 4 m but with probability p L =1− p S Ӷ p S ͒.The latter component is due to atmospheric turbulent-uplifting events, which have been up to now neglected in the prediction of invasion speeds.Both components are clearly distinguished in, e.g., Fig. 2 in Ref. 7 and they are due to independent events ͑because a seed can reach any final location either after being uplifted or not, but both possibilities are obviously incompatible͒.Therefore, the dispersion kernel can be written in the additive form ͑⌬ x , ⌬ y ͒ = p S S ͑⌬ x , ⌬ y ͒ + p L L ͑⌬ x , ⌬ y ͒.

III. CONTINUOUS-SPACE RANDOM WALK "CSRW… MODEL
For the evolution equation ͓Eq.͑2͔͒, assuming R 0 Ͼ 1 and that p͑x , y , t͒ has bounded support ͑i.e., vanishes outside a finite region͒ at some finite value of time, Weinberger presented a general approach to derive the front speed. 16,17But Weinberger's approach has been almost always applied to one dimension. 18,19Thus, it will be clarifying to deal briefly with the 2D case in this section.Also, we will use a much simpler approach than Weinberger's.
The speed of radially symmetric 2D front solutions can be found most easily by assuming that for t → ϱ, the front is approximately planar at scales much larger than that of individual dispersal events.We can then choose the x axis parallel to the local velocity of the front.Let c ϵ͉c x ͉ stand for this speed ͑c y = 0 in the local frame just introduced͒. 8,23We look for constant-shape solutions with the form p = p 0 exp ͓−͑x − ct͔͒ as x − ct → ϱ and, as usual, assume that the minimum speed is the one of the front 8,9 ͑we will check this assumption by means of numerical simulations in Fig. 1͒.Then Eq. ͑2͒ leads to the asymptotic ͑t → ϱ͒ speed of 2D fronts, where is the modified Bessel function of the first kind and order zero.We have related the dispersal probability per unit area ͑⌬͒ ͑i.e., into a rectangular area d⌬ x d⌬ y ͒ to that per unit length ͑⌬͒ ͑i.e., into a 2D ring of area 2⌬d⌬͒, because the kernel in Ref. 7 is ͑⌬͒.It is easy to see that We have assumed an isotropic dispersion kernel ͑i.e., that depends only on distance ⌬ ϵ ͱ ⌬ x 2 + ⌬ y 2 ͒ and obviously Previous papers on Reid's paradox 18,19 usually apply the corresponding 1D result instead of ͑3͒-͑5͒.In fact, in one dimension, Eq. ͑3͒ holds but Eqs.͑4͒ and ͑5͒ do not.Thus, speed c is different in one dimension than in two dimensions ͓indeed, in Sec.VI we show that the 2D speed is always slower than in one dimension for the same dispersal kernel ͑⌬͔͒.
The fact that 1D and 2D speeds are different was noted already in Ref. 22 for the evolution equation in ͑1͒ and in Ref. 10 for its hyperbolic and Fisher limits.However, here we will deal with the evolution equation in ͑2͒.
Since we are interested in Reid's paradox, which refers to forest range expansions that took place in two dimensions, For bimodal kernels it is found that computer simulations cannot yield accurate results within a reasonable computing time ͑Sec.IV͒.The DSRW model overcomes this limitation ͑Fig.2͒.
we will firstly deal with the 2D case.We shall compare speeds in two dimensions and one dimension only in Sec.VI.

IV. MOLECULAR DYNAMICS SIMULATIONS
We are not dealing with a differential but with an integrodifference equation in two dimensions, Eq. ͑2͒.Therefore, numerical simulations in this paper will not be based on finite-step approximations to derivatives but on molecular dynamics ͑or cellular automata in the continuous limit͒.Simulations in two dimensions are much more time costly than in one dimension, but a 2D model is necessary ͑as seen in the last paragraph of the previous section͒.We have thus performed simulations on a 2D grid, with nearest neighbors separated by a distance D. Initially p͑x , y ,0͒ =1 at ͑x , y͒ = ͑0,0͒, and 0 elsewhere.At each time step, we compute the new number density of trees p͑x , y , t + T͒ at all nodes of the 2D grid as follows.In agreement with Eq. ͑2͒, we first compute the seed production R 0 p͑x , y , t͒ at every node 24 and then redistribute this value among all grid nodes using the kernel ͑⌬͒.We have performed such 2D simulations for values of R 0 and T typical of the yellow poplar (Liriodendron tulipifera). 25We consider this species because its long-and short-distance kernel components ͓ S ͑⌬͒ and L ͑⌬͒, respec-tively͔ were determined in Ref. 7. Consider first a very simple, short-distance unimodal kernel S ͑⌬͒ such that it is approximately constant for dispersal distances ⌬Ͻ15 m and zero for ⌬Ͼ15 m. 7 Using a 2D grid with nearest neighbors separated by a distance D = 1 m, the simulations agree with the CSRW, as shown in Fig. 1. 26 This shows ͑i͒ the validity of the minimum-speed conjecture 27 and ͑ii͒ the need to take Eq.͑6͒ into account in the simulations. 28But S ͑⌬͒ is a unimodal kernel, whereas we are interested in bimodal ones ͑see Sec.I͒.Then the simulations above are not useful, as we have found that the time required is prohibitively long. 29herefore, molecular dynamics simulations are not practical to test the 2D analytical result in ͑3͒ for bimodal kernels.Is there some other way to test whether Eq. ͑3͒ holds or not for bimodal kernels?In the next section we present a fast, efficient approach that we have found very useful for this purpose.

V. DISCRETE-SPACE RANDOM WALK "DSRW… MODEL
This model is not exact, but it is necessary to check the CSRW model for bimodal kernels.The DSRW is closely analogous to the numerical simulations, in the following sense.Both in the DSRW and the simulations, we replace 2D continuous space by a grid of points ͑nodes͒ with nearest neighbors separated by a distance D along the x and y axes.The nodes are the only points available for seeds and trees.First consider the very simple, highly idealized case in which any tree disperses seeds only to its eight nearest neighbors on the grid.Obviously, these eight final dispersal nodes lie on a square with side 2D and center at the parent tree, as follows.The closest four nodes are a distance ±D away along the x or y direction, and the next four are at distance ±D away along both directions, i.e., on the vertices of the square ͑at distance where the first four terms correspond to horizontal and vertical "jumps," whereas the last four terms are due to diagonal jumps, and the jump probabilities are, from Eq. ͑6͒, For the simple case of Eq. ͑8͒, n = 2 and the only possible dispersal distances are ⌬ 1 = D and ⌬ 2 = D ͱ 2.
We apply the approach already used above ͓Eq.͑3͔͒ ͑Refs.8, 23, and 27͒ and obtain, from the DSRW equation ͓Eq.͑8͔͒, Note that Eq. ͑10͒ is a very simple approximation ͑DSRW͒ but is completely analogous to the exact ͑CSRW͒ speed in ͑3͒.We have found ͓e.g., for S ͑⌬͒ in Sec.IV͔ that this extremely simple DSRW yields a speed in ͑10͒ which disagrees with that from the CSRW.Thus, we next consider dispersal to nodes not on a single but on many squares ͑j =1,2,3, ...͒ centered at each parent tree.A square with side 2jD will obviously have 8j nodes, namely, four at dis-tance jD, four at distance jD ͱ 2, and also ͑except for the simple case j = 1 above͒ eight nodes at distance ͱ ͑jD͒ 2 + ͑iD͒ 2 for i =1,2, ... , j − 1.Finally, in order to use the measured kernels 7 we need to restrict dispersal to a maximum distance in whatever direction, r max .Then it is not difficult to write the analog to Eq. ͑8͒ for a bimodal kernel ͓i.e., p L L ͑r͒ + p S S ͑r͔͒ and see that the speed in ͑10͒ is generalized into

͑11͒
where N L = r max L / D L , N S = r max S / D S , and the terms with ͱ ͑jD S ͒ 2 + ͑iD S ͒ 2 arise from jumps in directions different from 0°, ±45°, 180°, and ±90°.The probabilities are related by Eq. ͑9͒ to the corresponding dispersion kernel, for example, For the yellow poplar ͑Liriodendron tulipifera͒, the LDD component of the kernel derived ͑and favorably compared to observations͒ in Refs.

· ͑13͒
so r max L =10 4 m, whereas, as mentioned in Sec.IV, its shortdistance component S ͑⌬͒ can be taken as approximately constant for ⌬Ͻ15 m and zero for ⌬Ͼ15 m, so r max S =15 m. Figure 2 presents the results for the bimodal kernel p L L ͑r͒ + p S S ͑r͒, where p L = 0.002 02 and p S =1− p L are the probabilities of long-distance and short-distance dispersals ͑obtained from Fig. 2 in Ref. 7͒.The results for the unimodal kernels L ͑r͒ and S ͑r͒ are also presented for comparison.The relevance of these results for Reid's paradox will be discussed in Sec.VII.

VI. COMPARISON OF INVASION SPEEDS IN TWO DIMENSIONS AND ONE DIMENSION
The aim of this paper has been to analyze front speeds in 2D space.But many papers have applied 1D results instead, see, e.g., Refs.18, 19, 30, and 31 and references therein.Thus, a comparison between two dimensions and one dimension is appropriate.
In one dimension, clearly the kernel per unit area ͑⌬͒ has no physical meaning.The 1D evolution equation equivalent to Eq. ͑2͒ is obviously where we have not used the notation ͑⌬͒ for the kernel in order to avoid confusion.Indeed, in two dimensions we have used Eq.͑7͒ for the normalization of the kernel per unit length , In contrast, in the 1D Eq. ͑14͒ the kernel per unit length ˜is defined both for positive and negative jumps.Therefore, the lower integration limit must now be −ϱ instead of 0, i.e., Clearly, ⌬ x can be either positive or negative, whereas ⌬ ജ 0 ͑they are related as ⌬ ϵ͉⌬ x ͉ in one dimension and ⌬ ϵ ͱ ⌬ x 2 + ⌬ y 2 in two dimensions͒.In order to compare to two dimensions, note that for an isotropic kernel, Eq. ͑16͒ implies that If we introduce the following 1D kernel ͑which we define only for ⌬ ϵ͉⌬ x ͉ ജ 0͒: which agrees with the corresponding normalization equation in two dimensions, see Eq. ͑15͒.This will make it possible to perform a meaningful comparison between one dimension and two dimensions.Now that we have carefully defined the kernel per unit length ͑⌬͒ both in one dimension and two dimensions, we can ask the following question.For a given kernel ͑⌬͒, is the front in one dimension faster or slower than in two dimensions?Up to now, this question seems to have remained unanswered for arbitrary kernels ͑⌬͒.Here, we will provide a mathematical proof that the front is always slower in two dimensions than in one dimension for the same kernel ͑⌬͒.
For an isotropic 1D kernel, i.e., ˜͑−⌬ x ͒ = ˜͑⌬ x ͒, Eq. ͑3͒ holds for the front speed c but ˆ͑͒ is not given by the 2D Eq. ͑4͒.Instead, in one dimension This result was first derived by Weinberger 16,17 and has been widely applied 18,19,30,31 ͑it can be also easily obtained from the approach in Sec.III above͒.The 1D speed can thus be written in terms of ͑⌬͒, from Eqs. ͑20͒ and ͑18͒, The 2D speed for an isotropic kernel ͑⌬͒ is, from Eqs. ͑3͒-͑5͒, where we have simply applied that cos͑␣ + ͒ = −cos ␣.Alternatively, Eq. ͑22͒ can be also obtained from Eq. ͑24.102͒ in Ref. 32.Finally, we write the 2D speed in ͑22͒ as where we have simply used the definition cosh x ϵ 1 2 ͓exp͑x͒ + exp͑−x͔͒.
We now rewrite the 1D speed in ͑21͒ in a form similar to the 2D speed in ͑22͒.Obviously,

͑24͒
Comparing Eqs.͑23͒ and ͑24͒ it is obvious that, for any symmetric kernel per unit length ͑⌬͒ ͑defined for ⌬ജ0͒, simply because cos ഛ 1 and cosh x ϵ 1 2 ͓exp͑x͒ + exp͑−x͔͒ increases with increasing values of x ജ 0 ͑in other words, its derivative is sinh x ϵ 1 2 ͓exp͑x͒ − exp͑−x͔͒ ജ 0 if x ജ 0͒.Therefore, we have shown that all symmetric kernels ͑⌬͒ give a speed which is lower in two dimensions than in one dimension.This result is not surprising, according to the following two arguments.

͑i͒
We may understand intuitively the result in ͑25͒ as follows.Consider, for example, a delta kernel, i.e., such that all individuals move the same distance.In one dimension, half of them will travel in each direction.In two dimensions, fewer individuals can reach a given point than in one dimension because there are not only two final dispersal points but a whole a circle.So, we could well expect the front to be slower in two dimensions than in one dimension.Moreover, any kernel can be regarded as a sum of delta kernels.
From this perspective, it is not surprising that fronts are slower in two dimensions than in one dimension.͑ii͒ For the special cases of Fisher and hyperbolic reaction-diffusion evolution equations, fronts are also slower in two dimensions than in one dimension. 10owever, it must be remembered that ͑i͒ such approaches are diffusive limits to Eq. ͑1͒, which does not take proper care of the time order of events ͓in contrast to Eq. ͑2͒ or ͑14͒, considered here͔, and ͑ii͒ diffusive limits are well known to break down for long-distance dispersal. 18l of the results in the present paper, including Eq. ͑25͒, have been obtained for the time-ordered evolution equations ͑2͒ in two dimensions and ͑14͒ in one dimension.Thus, they may not hold for other evolution equations ͓such as ͑1͔͒.
In this section we have shown that, for time-ordered evolution dynamics and, an arbitrary symmetric kernel per unit length ͑⌬͒, the front is always slower in two dimensions than in one dimension. 33or the short-distance kernel S ͑r͒, we include the predicted speed in ͑21͒ in one dimensional as a dotted curve in Fig. 1.It is seen that the 1D speed is faster than the 2D speed, as it should.If the 1D approach is applied, the error ͑difference relative to the 2D prediction͒ is important and increases rapidly for decreasing values of R 0 ͑for example, the error is 22% for R 0 =6 yr −1 and 33% for R 0 = 2.2 yr −1 ͒.This could have been expected since the lower the net reproductive rate, the more important is dispersal relative to reproduction.The differences between the 2D and 1D results show the importance of using 2D formulas to predict invasion speeds.Similarly, for the bimodal kernel p L L ͑r͒ + p S S ͑r͒, Fig. 3 shows the predicted speed in ͑21͒ in one dimension as a dotted curve.Again, it is faster than in two dimensions, as it should, and the differences are important, so it is not reasonable to use 1D models for invasions that take place in two dimensions ͑e.g., postglacial forest range expansions͒.

VII. DISCUSSION
Here, molecular dynamics simulations have been useful to ͑i͒ make us realize that ͑r͒ must be related to the measured kernel ͑r͒ ͑Ref.28͒ and ͑ii͒ check that the minimum speed in ͑11͒ is that of the front. 27Thus, it is not reasonable to rely on a single analytic approach only ͑e.g., the CSRW͒.Indeed, in reaction-diffusion problems, simulations are almost always used to check the analytic results. 8,9However, for bimodal kernels the 2D simulations do not yield accurate enough results within a viable computing time ͑Sec.IV͒.Thus, in addition to the 2D CSRW model, we have also presented a discrete analytical model ͑DSRW͒ and applied it ͑instead of the simulations͒ to check the CSRW results ͑Figs. 2 and 3͒.
Both the 2D DSRW and 2D CSRW models show conclusively that the front speeds for the bimodal kernels in Ref. 7 ͓i.e., p L L ͑r͒ + p S S ͑r͔͒ are about 10 2 -10 3 m / yr, i.e., two orders of magnitude faster than those for the unimodal, short-range component S ͑r͒ ͑Fig.2͒.Speeds of 10 2 -10 3 m / yr are, in fact, those required to solve Reid's paradox. 18hort-distance kernels S ͑r͒ have been measured experimentally many times.But bimodal kernels with a LDD component L ͑r͒ were derived by a mechanistic ͑or physical͒ model in Ref. 7 which motivated the work here reported.
We conclude that Reid's paradox of rapid tree migration can be solved ͑as far as the order of magnitude is concerned͒ by taking into account the bimodal dispersal kernels derived and favorably compared to data in Refs.7, 20, and 21.This is the main biophysical result of this work.
On the analytical side, our main contributions here are the DSRW model, taking into account two dimensions in the CSRW, and the comparison between two dimensions and one dimension ͑Sec.VI͒.But the important natural phenomenon tackled is Reid's paradox, which motivated all of our analytical work and 2D simulations.Of course, we do not claim to have solved Reid's paradox.In future work, it would be interesting to apply the 2D CSRW and DSRW approaches presented here by comparing the predictions and observations for a list of tree species such that ͑i͒ their range expansion speeds can be measured from the fossil record and ͑ii͒ dispersal kernels and net reproductive rates can be estimated using the methods in Ref. 7 and experimental measurements.Our results here do show the procedures ͑2D CSRWs and 2D DSRWs͒ useful to perform such an analysis ͑as well as the limitations of using numerical simulations for this purpose͒.Section VI contains the proof that, for an arbitrary kernel ͑⌬͒, the front is always slower in two dimensions than in one dimension under the time-ordered evolution equation ͑2͒ ͑in two dimensions͒ or Eq.͑14͒ ͑in one dimension͒.
It future work, it would be nice to apply the physical models presented ͑2D CSRWs and 2D DSRWs͒ to extended 2D evolution equations that take into account age structure, interspecific interactions, etc.It would be also interesting to generalize the models described here to a variety of applied biophysics systems, e.g., tumor growth models 1 with several nutrients ͑each with a different diffusion coefficient or dispersal kernel͒, several-species dispersion in microchannel flows, 2 etc.

FIG. 1 .
FIG. 1. Front speed in two dimensions vs net reproductive rate ͑Ref.25͒ for a unimodal short-distance kernel S ͑⌬͒.Stars: 2D computer simulations.Full curve: analytical 2D CSRWs, Eqs.͑3͒-͑5͒.There is good agreement.The 1D speed for the same kernel is included for comparison ͑dotted curve͒.For bimodal kernels it is found that computer simulations cannot yield accurate results within a reasonable computing time ͑Sec.IV͒.The DSRW model overcomes this limitation ͑Fig.2͒.

FIG. 2 .
FIG. 2. Front speeds in two dimensions vs net reproductive rate ͑Ref.25͒.Curves: CSRWs, Eqs.͑3͒-͑5͒.Symbols: DSRWs, Eq. ͑11͒, using the values of D L and/or D S in the legend ͑in meters͒ and the corresponding kernel͑s͒.The bimodal kernel for the yellow poplar, from Ref. 7, leads to the middle curve.It thus predicts speeds of about 10 2 -10 3 m/generation.In contrast, the short-range unimodal kernel ͑lower curve and stars, the same as in Fig. 1͒ predicts front speeds several orders of magnitude lower.Thus, the bimodal kernels derived and favorably compared to data in Refs.7, 20, and 21 may solve Reid's paradox.Note from the upper curve ͑100% of seeds with LDD, L ͒ that the same order of magnitude ͑10 2 m/yr͒ is obtained for only 0.2% of seeds with LDD ͑middle curve͒.

FIG. 3 .
FIG.3.Comparison of speeds in two dimensions and one dimension for the bimodal kernel.The full line and symbols correspond to two-dimensions ͑the same as in Fig.2͒, whereas the dotted curve corresponds to one dimension.The 1D speed is about 20% faster than the 2D speed.