Diffusion processes are commonly used for modelling natural space-time processes in climate sciences, ecology, epidemiology and population dynamics, to name a few. For example, the diffusion equation, with , one of the most elementary partial differential equations (PDEs) known as the heat equation, describes the macroscopic behaviour of many micro-particles resulting from their random movements. How PDEs such as the heat equation, but also more complex ones, interact with spatial statistics is the focus of this contribution.
PDEs and statistics have long been connected through the usual question of parameter estimation in mechanistic models (see e.g., Wikle, 2003; Soubeyrand and Roques, 2014, and the references therein). About ten years ago, Lindgren et al. (2011) revisited the explicit link between stochastic partial differential equations (SPDEs) and spatial statistics already established in Whittle (1954, 1963). They introduced a new paradigm for the statistical analysis of spatial data that consists in modelling a spatial variable as the solution of a certain class of SPDE for which the representation on a finite grid presents interesting Markov properties, making it computable for very large grids and capable of handling very large datasets. This approach was based on the observation made in Whittle (1954, 1963) that, on , the solution of the SPDE
is a Gaussian random field whose covariance function is proven to be a Matérn covariance. In (1) , is the usual Laplacian,
is a Gaussian (spatial) White Noise process andis the unknown stationary random field. For , the pseudo-differential operator in (1) was constructed on an ad hoc basis to get a Matérn covariance. As Whittle (1954) himself indicated “it is difficult to visualize a physical mechanism which would lead to (1)”. However, in contrast to the previous statistically oriented constructions, the SPDE approach should also lead to a physically grounded construction, for which the parameters would carry traditional physical interpretation such as diffusivity, reaction and transport. The SPDE approach also allows the construction of models with interesting properties, such as non-symmetry and space-time non-separability (Brown et al., 2000). Non-stationarity can also easily be accounted for by letting the parameters in (1) be space-dependent, see e.g. Lindgren et al. (2011) and Fuglstad et al. (2015). Space-time generalisation is not straightforward. In particular, the space-time covariance is not always accessible (but the spectral density is) and the definition of the space-time white noise requires some care. Recent advances have been made in Bakka et al. (2020) and Carrizo Vergara et al. (2021), where space-time models admitting Matérn covariance functions in space have been proposed.
Carrizo Vergara et al. (2021) proposed a rigorous mathematical framework, based on the theory of Generalized Random Fields (Itô, 1954; Gelfand and Shilov, 1964), to construct spatio-temporal random fields arising from a very broad class of linear SPDEs. They considered the general class of linear SPDEs of the form defined through operators of the form
where, in the real case, is a continuous, symmetric function bounded by a polynomial,
the Fourier transform onand the inverse Fourier transform. The function in (2) is referred to as the symbol function of the operator. If, in addition, is bounded from below by the inverse of a strictly positive polynomial, there exists a unique Generalized Random Field solution to the equation where ‘2nd. o.’ means that the solution has to be understood in the second order sense (see Carrizo Vergara et al., 2021, for more details). In this case, the spectral measures of and are connected by:
which leads to the covariance . Thus, the relationship between the covariances of and is
When reduces to a white noise , its spectral measure is , and the covariance reduces to the Dirac . Hence, the covariance can be seen as a Green’s Function of the operator . This construction offers the possibility to build and characterize models far beyond the Matérn family which is currently the covariance model considered within most SPDE implementations.
In this contribution, our aim is to further explore the intimate interplay between partial differential equations, their stochastic extensions, namely the SPDEs, and spatial statistics. Most of the literature about PDEs focused on elliptic and parabolic operators involving a Laplacian or its anisotropic extension. Although the fractional Laplacian has recently generated much attention in the PDE community, its microscopic interpretation is much less classical than that of the standard Laplacian. As a preliminary to this study, we propose in Section 2 an introduction to the connection between Lévy flights and PDEs involving the fractional Laplacian operator. The corresponding Fokker-Planck PDEs will serve as a basis to revisit in Section 3 the original SPDE framework associated with Eq. (1), but with a mechanistic perspective. We propose new generalisations by considering a general form of SPDE with terms accounting for dispersal, drift and reaction. Regarding the dispersal term, we first show that there is an equivalent integral representation offering a mechanistic interpretation of the operator in the left-hand-side of Eq. (1). As a side result, we show the profound difference between the fractional Laplacian , which is associated with a fat-tailed dispersal kernel (for ) and therefore describes long-distance dependencies, and the damped fractional Laplacian , which is associated with a thin-tailed kernel, corresponding to short-distance dependencies. We then introduce the fractional Laplacian operator with a linear reaction term and we also consider kernel-based dispersal operators. In each case, we derive the associated covariance function and we compare its properties with those of the Matérn family. We then illustrate SPDEs with non-stationary drift terms representing external spatially and temporally varying forces. Finally, we introduce nonlinear bistable reaction terms ; we discuss their physical meaning in a deterministic framework and the possible applications of the corresponding SPDE, the stochastic Allen-Cahn equation. Lastly, in Section 4, we revisit the classical use of statistical models for parameter estimation in PDE models, with which we started this introduction. In this original approach, we combine the Fokker-Planck PDEs of Section 2 with inhomogeneous point processes. Section 5 concludes with some elements of discussion, with a particular focus on the form of the diffusion operator.
In the following, will represent a random field defined over space and time, with and . Partial derivatives along time will be denoted and is the Laplacian. The Fourier transform is defined, as follows: . As already mentioned in the first paragraph, we will often adopt a microscopic point of view, in which solutions of PDEs or SPDEs describe the density of particles subject to diffusion equations.
2 Lévy flights and PDEs with fractional Laplacian
Since Einstein’s theory of Brownian motion (Einstein, 1905), the parabolic PDE has a well-established physical interpretation. The solution can indeed be seen as the density of particles following independent Brownian motions. More generally, Fokker-Planck equations act as a natural bridge between parabolic PDEs and Itô diffusion stochastic differential equations (SDEs) (e.g., Gardiner, 2009). These SDEs extend Brownian motion by introducing a drift (or similarly a transport) coefficient and spatio-temporal heterogeneities in the movement of particles. The solution of parabolic-like equations but with a fractional Laplacian (with ) instead of a standard Laplacian () can also be interpreted as a density of particles. But in this case the particles follow independent Lévy flights. As this physical interpretation is less standard, we propose below a short introduction to the connection between Lévy flights and Fokker-Planck PDEs involving a fractional Laplacian.
We consider a particle whose position follows a dimensional stable Lévy process, for some (see e.g., Cont and Tankov, 2003; Applebaum, 2009, for precise definitions). This can be described by a Langevin-like equation (e.g. Schertzer et al., 2001; Applebaum, 2009):
where is an stable rotation-invariant -dimensional Lévy process with scale parameter and is a forward increment in time of this process. As mentioned above, when , is a standard Brownian motion and the equation (3) is a standard Itô diffusion SDE (for a discussion of the differences between Itô and Stratonovitch SDEs, the reader can refer to Gardiner, 2009; Horsthemke and Lefever, 2006)
. The vector fieldis the drift coefficient. It is related to the existence of some external force acting on the particles, sometimes also called bias, e.g. due to wind or attraction/repulsion. For simplicity, we assume an isotropic and time-independent random driving force (but see Remark 1 below), which means that the term (the diffusion matrix when ) satisfies where is a scalar function and
is the identity matrix of size. Note that, in general, the coefficients and may depend on the current position . For example, the particle motion could be slower in some regions of space, which would be modelled with a smaller value of in these regions. Note that the existence and uniqueness of the solution of (3) are guaranteed under some technical growth conditions on the coefficients (see e.g., Arnold, 1974; Schertzer et al., 2001; Øksendal, 2003).
The Fokker-Planck equation describes the dynamics of the transition probability density. We denote by this transition probability. In the standard diffusion case , the dynamics of are described by the parabolic PDE (see e.g. Gardiner, 2009):
where the divergence operator ‘div’ acts on the vector fields , When the movement of the particles is not driven by a Gaussian noise (i.e. when ), the Fokker-Planck equation involves a fractional Laplacian (Schertzer et al., 2001):
see Section 3 for more details about this operator; see also Kwaśnicki (2017) for alternative definitions of this operator. Notice that, with and (Brownian diffusion), the equation (5) is equivalent to (4).
In the general anisotropic case, when and is not necessarily of the form , we obtain an equation similar to (4), but the term must be replaced by with the entries of the matrix . When does not depend on , we note that this operator reduces to . Otherwise, it can be written as a sum of a ‘Fickian’ diffusion term (see also the discussion in Section 5) and a drift term: , with the vector with coordinates .
The transition probability density involves an initial condition , where is a Dirac measure at . Assuming independent particles, with common initial distribution , it is then straightforward to check that the expected density of particles satisfies the PDE (5), with the initial condition
Assume further that the particles have a life expectancy for some and disappear at an exponential rate. Consider any given particle. We denote by its “death” time (or removal time: we not only focus on living organisms). With a slight abuse of language, we define as the probability that the particle is at the position and still dispersing at time , that is, for any measurable set ,
Since dispersal and death are defined as two independent processes, . Therefore, and consequently the expected particle density satisfy the following equation:
This construction highlights two interesting connections between PDEs and stochastic processes. On the one hand, it illustrates that some PDEs with fractional Laplacian arise quite naturally as the Fokker-Planck equation of SDEs with Lévy flights. On the other hand, the PDE in (6) suggests interesting generalisations of the SPDE in (1) which will be the subject of the next section.
3 SPDEs and random fields: a mechanistic point of view
Inspired by (6), we consider a quite general class of time-dependent SPDE models that rely on mechanistic assumptions concerning the underlying variable. These models take the following general form:
The random field is thus governed by a dispersal term , a drift term , a reaction term describing the local feedback of the random field on itself and a spatio-temporal white noise defined as the derivative of a Wiener process, and thus verifying . The square bracket symbol is used to underline that the operator may depend on the whole function , and not only on its value at . The well-posedness of (7) depends on its precise form and on the space dimension . When the reaction term is linear and with coefficients that are constant in , the existence and uniqueness of a stationary distribution-valued solution generally follow from the results in Carrizo Vergara et al. (2021), at least in the examples that we treat below. With nonlinear reaction terms (nonlinear with respect to ), the theory is also well-established in the case (e.g. Da Prato and Zabczyk, 2014). The picture is less clear when is nonlinear and : although such equations are regularly used in applied sciences, mathematical results tend to show that additive white noise leads to ill-posed equations (Hairer et al., 2012; Ryser et al., 2012).
In this equation, there is a diffusive dispersal term and an absorption term (the reaction) which forces the solution to remain close to , see Section 2 for a microscopic interpretation. When the noise does not depend on time (spatial white noise ), the SPDE (1) in Lindgren et al. (2011) is a time-independent solution of () for . In the general case the SPDE (1) involves the “damped” fractional operator , which does not exactly fit in (7).
We will discuss below three extensions of the “standard” SPDE-GMRF that go beyond the framework () set in Lindgren et al. (2011) and their possible applications in environmental and ecological sciences. These extensions either deal with the dispersal term (Section 3.1), the drift term (Section 3.2) or the reaction term (Section 3.3).
3.1 Dispersal operators
In equations of the form (7), nonlocal dispersal operators
can describe short or long-distance dispersal, depending on the precise shape of the kernel (Klein et al., 2006; Garnier, 2011). When imposing that the integral of is equal to 1, the quantity can be interpreted as the probability density that a dispersing particle initially located at moves to the position (see e.g., Roques, 2013), and the coefficient is the rate of dispersal. The frequency of long-distance dispersal events is governed by the tail of the kernel as : thin-tailed kernels decay at least exponentially and fat-tailed kernels decay more slowly than any exponential function.
Equations of the form (7) with dispersal terms of the form (8) and regular kernels of mass , null drift term and fit in the class of SPDEs considered in Carrizo Vergara et al. (2021), as , with symbol function , which is an Hermitian-symmetric function (see Theorem 1 in Carrizo Vergara et al., 2021). We also consider singular kernels for which, even if is not properly defined, the right-hand side in (8) is still well-defined, for instance when . We compare here, from a mechanistic viewpoint, the physical interpretation of these equations
compared to evolving Matérn equations
The main findings of this section are summarized in Table 1.
3.1.1 Fractional Laplacian
We begin by recalling that the fractional Laplacian has a well-known pointwise representation when (e.g., Kwaśnicki, 2017): for all (),
for some constant . Due to a singularity around , the above integral has to be understood in the Cauchy principal value sense when . Equation (11) shows that the fractional Laplacian is a nonlocal operator, and formally corresponds to a dispersal term of the form (8), although the kernel is not integrable around . The kernel belongs here to the family of (very) fat-tailed kernels, which is consistent with the microscopic interpretation of this operator as seen in Section 2, since Lévy flights (3) are continuous time random walks, with possible large jumps if .
3.1.2 Fractional damped Laplacian
We recall that the solution to in dimension is a Gaussian random field with Matérn covariance of parameter :
where is the modified Bessel function of the second kind with index . Notice that in the limit case , is a generalized Matérn covariance which is well defined for .
In order to find kernel associated to , we use a semigroup formula (e.g., Stinga, 2019) (for and , ):
with the gamma function. The quantity is the semigroup generated by acting on , i.e., the solution of the Cauchy problem:
More explicitly, , with the heat kernel. Replacing into (13), we have:
Then, we note that:
Finally, using Fubini’s theorem and assuming the convergence of the double integral in (14), we formally get the following pointwise expression for the operator , in the Cauchy principal value sense:
This formula leads to a mechanistic interpretation of the operator : it contains an absorption term (the reaction term ) which becomes stronger as is increased, and a nonlocal dispersal term of the form (8). The kernel behaves as for large : it is therefore a thin-tailed kernel, and its tail becomes thinner as is increased. Contrarily to the fractional Laplacian (with ), and similarly to the standard Laplacian, the operator therefore describes short-distance dispersal events. If we come back to the definition of the Matérn covariance (12), we note that has an exponential decay (like ) as , corresponding to short-distance dependencies.
3.1.3 Fractional Laplacian with linear reaction
The above analysis shows that the evolving Matérn equation (10) can be written in the form (9), with , and replaced by . This emphasises the interest of using equations of the general form (9) with kernel operators of the form (8) as they should lead to additional flexibility, compared to evolving Matérn equations. We consider here the Eq. (9) with , which amounts to considering the fractional Laplacian operator with linear reaction term:
in the time-dependent and time-independent cases (when does not depend on , and with a spatial white noise ), respectively. In the time-independent case, some explicit expressions for the stationary covariance function can be obtained. First, with , and in dimension , we apply a Fourier transform to obtain the stationary covariance function:
with . We note that is the so-called “first auxiliary function” associated with trigonometric integrals, which can also be written (Chapter 5 in Abramowitz and Stegun, 1964), with and the Cosine and Sine integrals. Using the asymptotic expansions and , we get as . With the same parameter values (, ), recall that the damped Laplacian leads to the 1/2-Matérn covariance function, which decays exponentially for large .
In dimension , still with , we use the relationship between the Fourier transform of a radially symmetric function and the Hankel transform to compute the stationary covariance:
with the Bessel function of the first kind. The last integral in the above equality is a tabulated integral, and corresponds to the Hankel transform of order 0 of : from (case 126.96.36.199 in Prudnikov et al., 1986) we know that
with the Bessel function of the second kind and the Struve function of order . Differentiating this expression with respect to and using (18), we get:
Here, as . Note that, with and , we have In this case, the damped Laplacian does not lead to a Matérn covariance, but to a generalized covariance function proportional to as (Carrizo Vergara et al., 2021; Allard et al., 2021).
In both cases and , we have thus observed that the covariance function associated with the fractional Laplacian with absorption has an algebraic decay as . As expected from the mechanistic interpretation of the operator, the obtained covariance function therefore describes long-distance dependencies. See Fig. 1 for an illustration of the behaviour of the new covariance function (19), compared to the generalized covariance given by the damped fractional Laplacian in the case and .
3.1.4 Convolution kernels
We now move away from pseudo-differential operators and consider the equation (9) with regular convolution kernels with mass , for which is well-defined in (8). The spectral measures are (Theorem 1 in Carrizo Vergara et al., 2021):
in the time-dependent and time-independent cases, respectively. As as , these measures are not finite, see Remark 1 in (Carrizo Vergara et al., 2021). In the time-independent case, with an exponential kernel , we compute the Fourier transform of thanks to the Hankel transform:
where the last equality follows from (Prudnikov et al., 1986, case 188.8.131.52 with and ). In the case ,
Since and (for ), applying a Fourier transform, we obtain an expression for the stationary covariance function:
for some positive constants , . In this case, we observe that the covariance function is comparable to a 3/2-Matérn covariance with an additional Dirac mass at , referred to as the “nugget effect” in geostatistics.
Let us now consider the more general case where the Fourier transform of is , with
. This corresponds to an exponential kernel with an odd value of the dimension, see Eq. (20) or more generally to a dispersal kernel given by a Matérn covariance in with parameter . Then,
The spectral measure associated to a general Matérn dispersal kernel is thus the sum of a Lebesgue measure and two spectral densities which are the inverse of positive polynomials of . Thanks to the Rozanov Theorem (Rozanov, 1977), it is known that the Gaussian random fields on characterised by such spectral densities are Gaussian Markov random fields, and thus with thin-tailed covariances.
As the kernel were thin-tailed here, the exponential decay of the covariance functions at infinity is not surprising. We expect a different behavior for fat-tailed kernels (see Table 1 in Klein et al., 2006), leading to long-range dependencies as those observed with the singular kernel corresponding to the fractional Laplacian.