Speed of travelling fronts: Two-dimensional and ballistic dispersal probability distributions

The speed of traveling fronts for a two-dimensional model of a delayed reaction-dispersal process is derived analytically and from simulations of molecular dynamics. We show that the one-dimensional (1D) and two-dimensional (2D) versions of a given kernel do not yield always the same speed. It is also shown that the speeds of time-delayed fronts may be higher than those predicted by the corresponding non-delayed models. This result is shown for systems with peaked dispersal kernels which lead to ballistic transport.

Introduction. -Wavefronts are spatial profiles of concentration, temperature, etc., which connect two equilibrium states and travel without changing their shapes. Fedotov [1] has considered a kernel in which all particles jump the same distance. This highly idealized case yields very interesting comparisons between alternative evolution equations. Kot and coworkers [2] have derived important results (e.g., accelerating fronts and approximate profiles) for integrodifference evolution equations. In biological invasions, such as animal migrations and the geographic spread of epidemics, it is usual to describe these processes by employing 1D reaction-diffusion models [3]. However, most real processes of biological invasions take place in 2D and, as we show in this work, the modelization of a 2D invasion cannot be analyzed with 1D models if one has to deal with dispersal probability distribution functions (dispersal kernels). We establish how the non-trivial correspondence between a 2D and a 1D model must be performed.
In this work we will start from a 2D continuous-time integral equation for the reactiontransport of particles which includes memory (or delay) effects described by a characteristic time T . Our model describes a more realistic situation than the classical 1D models because particles may jump in 2D according to a 2D dispersal kernel and react at the same time.
This is a potentially interesting model which should be taken into account to model biological invasions rather than the 1D models vastly employed in the literature. We will analyze the speed of fronts which emerge from a localized initial condition as we already did for 1D models [4]. Detailed results for 2D Gaussian and Laplace (both are of ecological interest) distribution kernels will be presented. We will show that Gaussian 1D distribution kernels yield the same speed as that for the 2D case. However, for Laplace distribution kernels the speed is different in 1D than in 2D. This is a result of a great ecological relevance. Numerical simulations for the molecular dynamics of the reaction-dispersal process for 1D and 2D jumps confirm all of our analytical predictions. Finally, we shall show that in time-delayed phenomena [5], ballistic-diffusive transport can lead to higher speeds than those predicted by the corresponding non-delayed models.
Reaction-dispersal model. -We derive the reaction-dispersal equation according to the continuous-time random walk (CTRW) framework. The quantity which defines the motion is the probability distribution Ψ(x, t) of the particle performing a jump of length |x| after waiting a time t at its starting point. If P (x, t) is the probability density of arriving at point x at time t and ρ(x, t) is the probability density of being at point x at time t, we have where φ(t) is the probability of remaining at least a time t on the point before proceeding with another jump. If ψ(t) = dxΨ(x, t) is defined as the waiting time distribution function, by the definition of φ(t) one has φ(t) = ∞ t dt ψ(t ). The term F [ρ(x, t)] stands for the probability density of the appearance of a new particle at point x at time t. As is usual in population dynamics, this term depends explicitly on the density ρ(x, t) in a non-linear form. The Fourier-Laplace transform of (1) yields whereF denotes the Fourier-Laplace of F [ρ(x, t)] and sφ(s) = 1 − ψ(s). If the jump length and waiting time are independent random variables, one finds the decoupled form Ψ( We assume that all the particles wait the same time T before performing a new jump (ψ(t) = δ(t − T )). Making use of the Taylor series and inverting by Fourier-Laplace eq. (2) becomes We make use, at this stage, of marginal stability analysis [6] in order to compute the selected speed of the front. Accordingly, the asymptotic speed v of the solution ρ(x, y, t) ∼ exp[i(k x x + k y y − ωt)] can be obtained from the dispersion relation ω = ω(k x , k y ) of the linear part of eq. (3) around the unstable state ρ = 0 (with ω, k x and k y complex in general) assuming Fisher-KPP kinetics (F [ρ(x, t)] F (0)ρ, where F (0) = (dF/dρ) ρ=0 ). The linear speed-selected v = ω * /k * satisfies k * = ik * i = i(k * ix , k * iy ), k * r = 0 and ω * = iω * i , ω * r = 0. The subscripts i and r indicate the imaginary and real parts, respectively. The value of ω * i must be obtained from the demand that the speed must have a minimum at ω * i . Then, the dispersal relation is where β ≡ F (0)T and We assume now that the dispersal process is isotropic and homogeneous, so that 2D dispersal kernels will depend explicitly on the radial jump length r only where r 2 = x 2 + y 2 . Then the integrals in (5) may be performed by using x = r cos θ, y = r sin θ to yield where k ≡ k 2 x + k 2 y , J 0 (kr ) is the Bessel function of the first kind and the normalization condition is then 2π ∞ 0 r dr ϕ(r ) = 1 and we have assumed radial symmetry in the second equality of (6). Note that eq. (6) is nothing but the Hankel transform of the 2D dispersal kernel. If one works with a 1D dispersal kernel, the dispersal relation is formally given also by eq. (4) but the integral transform for the kernel is basically a Fourier transform [4]. Therefore, in order to have the same speed of fronts for 1D and 2D models, the Fourier transform of the 1D kernel must be equal to the Hankel transform of the 2D model. The dispersion relation ω * i (k * i ) cannot be solved explicitly from eq. (4), so the speed must be found as v = min where k * i (ω * i ) is found form eq. (4). We will detail the calculations in the following section. Examples. -G a u s s i a n a n d L a p l a c e k e r n e l s. Gaussian and Laplace dispersal kernels have been widely used to fit the dispersal field data for D. pseudoobscura [7] and for the gypsy moth fungal pathogen [8]. Let us, in this section, illustrate how to derive the speed of the front for 2D Gaussian and Laplace kernels. We will show that, whereas for Gaussian kernels the speed derived in 1D and 2D is the same, for 1D and 2D Laplace kernels the speeds are different. This will be checked by means of numerical simulations in 1D and 2D. We start with Gaussian kernels. For a 1D Gaussian kernel, So that, the dispersal relation given by eq. (4) will be the same for 1D and 2D Gaussian kernels and will yield the same speed for the front. However, this does not occur between 1D and 2D Laplace kernels. For a 1D For this case, we observe that the dispersal relation (4) will be different for 1D and 2D Laplace kernel. We simulated the reaction-dispersal equation on a two-dimensional 2000 × 2000 matrix starting from a central origin where we assign the value 1, while the rest of the points were taken as 0. Every iteration consists of the following steps: for the dispersal, every non-zero point was considered as a new diffuser and its influence over its neighbors was determined by applying the kernel (i.e., eq. (3) without the second term in the right-hand side); at the same time, the logistic term F (ρ) = βρ(1 − ρ)/T up to third order (n = 3 in eq. (3)) was applied to all points to reproduce the reaction process (i.e., eq. (3) without the first term in the right-hand side). We observe in fig. 1 that slight deviations from the theoretical predictions appear as β/T grows; this is due to the discretization of the kernel necessary to perform the simulations in a discrete system. The discretization entails a loss of information, specially at the tail of the kernel. As the importance of the reaction process increases (it is, for high β/T ), the importance of the tail increases too; then, the discretization effects (plus restrictions concerning the running time of the simulations) make our simulations less accurate for high values of β/T . It is noteworthy to mention that the number of particles in our simulations has been taken large enough to avoid cut-off effects, shown before to involve a non-negligible decrease in the speed of reaction-diffusion wavefronts [9] regarding both discrete and continuous systems. Although these cut-off effects have not been considered here, we think that they represent a challenging field for further works.
The results obtained were travelling functions with a typical wavefront profile which spread in the radial direction. Due to the spatial homogeneity and the isotropy of our kernels, the level curves ρ = const are concentric circles in the plane XY with radius R(t) = vt for large times, where v is the radial speed of the front. To find the speed of these functions we measured the displacement of the mean height of the front per unit time and confirmed that, for sufficient long times (i.e., after a large number of iterations), it approached a constant. In fig. 1 we compare the results of these numerical simulations of the molecular dynamics of the process with the analytical predictions given by eq. (7). Note that for any value of β/T , the Laplace kernels always yield a higher speed than the Gaussian kernels. This is due to the shape of the tail of the kernel: the Gaussian kernel decays faster than the Laplace kernel and therefore leads to fewer jumps with very long lengths. Furthermore, it is seen that there is good agreement between the simulations and the analytical predictions from the model (3), and it is clear that 1D and 2D kernels do not always yield the same speed.
2D d o u b l e-d e l t a k e r n e l s. We will here derive a previously unknown, relevant and somehow surprising property of the new exact solution (4), which is at variance with the  predictions obtained for small T up to n = 2 in eq. (3). Consider first the still more restrictive, non-delayed approximation, namely is the diffusion coefficient. Equation (8) is nothing but the well-known Fisher-KPP value [3,10,11]. The superindex (1) means that this result can be also derived by taking up to n = 1 in (3). It is clear and well known [12] that the hyperbolic approach v (2) = 2 F (0)D(1+F (0)T /2) −1 (eq. (3) up to n = 2) always predicts a lower front speed than that given by the non-delayed approximation (8), v (2) < v (1) . In contrast, as we will now show, the exact equation (4) can lead to a higher front speed, v, than that predicted by the Fisher-KPP result (8), namely v (1) .
For our purposes now, it will be sufficient and much clearer to deal with a very simple case. Consider the 2D double-delta kernel where δ 2D αi (r ) is the 2D Dirac delta centered at r = α i , and p 1 and p 2 are positive numbers such that p 1 + p 2 = 1. The Fourier transform in eq. (9) can be easily found: where I 0 is the modified Bessel function of the first kind. Consider first p 1 = 0.1 and p 2 = 0.9 ( fig. 2a). Then, for the kernel in fig. 2a, T = 0.2 and F (0) = 1, the non-delayed or Fisher-KPP value (8) is v (1) = 21.2 and the exact, time-delayed result from eq. (4) is v = 19.7, so that v < v (1) . But in the case p 1 = 0.9 and p 2 = 0.1 ( fig. 2b), the exact result from eq. (4) is v = 7.7, whereas v (1) = 7.1, thus v > v (1) . This surprising result may seem counterintuitive at first sight. How may a delay lead to a higher speed? The answer is that the reason does not lie in the delay, but in the shape of the dispersion kernel. It is the combination of both the delay and the kernel shape effects which determines the speed. In order to see this, note again that Fisher's speed (8) depends on ϕ(x , y ) only through r 2 . For the kernel (9), r 2 = p 1 α 2 1 + p 2 α 2 2 . Therefore, a single 2D Dirac delta-function δ 2D γ (r ) will yield the same Fisher speed (8) (and also the same hyperbolic speed v (2) ) as the kernel (9), provided that γ is chosen such that r 2 is the same for the kernel (9) and for the kernel δ 2D γ (r ), i.e. γ 2 = p 1 α 2 1 + p 2 α 2 2 . As shown in fig. 2a, for the first example above, δ 2D γ (r ) (dashed curve) is close to the upper end of the range of the kernel (9). However, in the second example ( fig. 2b, dashed curve), δ 2D γ (r ) corresponds to a much lower dispersion distance than that of the particles to the right. These particles spread much further away than those described by the kernel δ 2D γ (r ). Thus, they may be expected intuitively to yield a higher front speed, which is the result obtained above (v > v (1) for the case in fig. 2b, whereas v < v (1) for fig. 2a). It is straightforward to check from eq. (4) that, for all of the single-peaked distributions in fig. 2 (dashed curves), one has v < v (1) . Thus we see that the result v > v (1) arises because of the double-peaked structure of the kernels in figs. 2b and c (full curves), in which a small part of the particles disperse to much higher distances than the rest. Therefore, it is the kernel shape that makes it possible for a time-delayed approach (eq. (4)) to yield a higher speed than the non-delayed one (eq. (8)), in spite of the fact that the delay time T tends to lower the speed. This may be also checked by noting that v (2) is time-delayed, but does not take into account the detailed shape of the kernel so it always yields v (2) < v (1) . The kernel shape is taken into account when the exact time-delayed approach (4) is used, because it includes infinite terms not only in the hyperbolic approach v (2) but also in the dispersaldistance expansion. The situation considered is somehow analogous to a similar one which arises in the kinetic theory of gases [13], radiative transfer theory and neutron transport [14], where one distinguishes between diffusive and ballistic behavior. For example, in the case of radiation inside an absorbing medium, diffusive transport corresponds to the photons emitted inside the medium, whereas the ballistic component is due to photons originating from the interface and reaching a point inside the medium without being absorbed (see fig. 2 in ref. [13]). Rather similarly, in the kernel presented in fig. 2c the decreasing component corresponds to the diffusive component, whereas the ballistic one is that on the right and consists of particles which make longer jumps.
Our new results can also be relevant in biophysical applications. For example, a basic unsolved problem is the modelling of the range expansion of species with several modes of dispersion. On the one hand, one has a local process of spread that can be modelled according to Fisher's equation (8) or its second-order extension v (2) (see, e.g., ref. [12]). But it is also observed experimentally that many species have a ballistic component of spread, mostly due to human-mediated dispersal which causes a small proportion of individuals to jump much higher distances, therefore yielding a higher front speed [15]. This corresponds precisely to our case, as depicted in fig. 2. We think that our theory can be very useful because no mathematical approach to these phenomena has been previously advanced.
Conclusions. -We have presented a general procedure to deal with ballistic-diffusive effects in front propagation. First of all, we have derived the speed of the front when 2D dispersal kernels are taken into account and we have stressed that in general 1D and 2D models do not yield the same speed even for isotropic kernels depending on the radial dispersal length ( fig. 1). Comparison with molecular dynamic simulations exhibits a very good agreement both for 1D and 2D dispersal kernels with our analytical predictions ( fig. 1). Secondly, we have come to the conclusion that the shape of the dispersion probability-distance distribution may lead to time-delayed fronts which travel at speeds higher than those corresponding to usual (non-delayed and hyperbolic) descriptions. This has been illustrated for several examples ( fig. 2), which have made it possible to understand the origin of this new effect.
Diffusive-ballistic time-delayed transport has been shown to arise naturally from the microscopic equations of kinetic theory, both for matter [13] and radiative [14] systems. This includes thermal conduction [13] in addition to diffusion. But, as explained at the end of the former section, biological invasions are a topic in which we think that our results can become very useful because the field dispersal data should be fitted to 2D kernels instead of 1D kernels as is common in the literature [7,8]. Likewise, another important field which can benefit from our approach is the study of extreme values in tree-like or network structures. More specifically, it has been shown recently that these systems can be analyzed by means of front solutions, exhibiting a wide range of potential applications to real problems (see [16]). * * * We thank D. Jou for pointing out ref. [13]