Finite-velocity diffusion on a comb

TL;DR

This paper analyzes a Cattaneo equation on a comb structure, deriving solutions for fractional diffusion that describe the transition from normal diffusion to subdiffusion due to the geometry.

Contribution

It provides a rigorous analysis of fractional diffusion equations on a comb, with solutions expressed via Fox H-functions, capturing the transition from normal to subdiffusive behavior.

Findings

Solutions in Fox H-function form for probability distribution

Mean square displacement expressed as an infinite series

Transition from normal diffusion to subdiffusion due to geometry

Abstract

A Cattaneo equation for a comb structure is considered. We present a rigorous analysis of the obtained fractional diffusion equation, and corresponding solutions for the probability distribution function are obtained in the form of the Fox $H$-function and its infinite series. The mean square displacement along the backbone is obtained as well in terms of the infinite series of the Fox $H$-function. The obtained solutions describe the transition from normal diffusion to subdiffusion, which results from the comb geometry.

Equations

$$\frac{\partial }{\partial t}u(x,t)+\tau \frac{{\partial }^2}{\partial t^2}u(x,t)=\mathcal{D}\frac{{\partial }^2}{\partial x^2}u(x,t),$$

$$⟨x^2(t)⟩\simeq t^{1/2}.$$

$$\frac{\partial }{\partial t}P(x,y,t)+\tau \frac{{\partial }^2}{\partial t^2}P(x,y,t)$$

$$P(x,y,t=0)=\delta (x)\delta (y),\frac{\partial }{\partial t}P(x,y,t=0)=0,$$

$$H_{p,q}^{m,n}(z)$$

$$=$$

$$\theta (s)=\frac{{\prod }_{j=1}^m\Gamma (b_j+B_{j}s){\prod }_{j=1}^n\Gamma (1-a_j-A_{j}s)}{{\prod }_{j=m+1}^q\Gamma (1-b_j-B_{j}s){\prod }_{j=n+1}^p\Gamma (a_j+A_{j}s)},$$

$$\tilde{\hat{P}}(k,s)=\frac{1+s\tau }{s(1+s\tau )+\mathcal{D}{k}^2}.$$

$$\hat{P}(x,s)$$

$$\hat{P}(x,s)=\frac{1}{2v}\sqrt{z}{e}^{-\rho \sqrt{z}}.$$

$$\hat{P}(x,s)={\mathcal{M}}^{-1}\left[\mathcal{M}[\hat{P}(x,s)](\xi )\right](s)$$

$$=\frac{1}{2v}{\mathcal{M}}^{-1}\left[{\int }_0^{\infty }{z}^{1/2+\xi -1}{e}^{-\rho \sqrt{z}}dz\right](s)$$

$$={\mathcal{M}}^{-1}\left[\frac{{\rho }^{-2\xi -1}}{v}\Gamma (2\xi +2)\right](s).$$

$$\hat{P}(x,s)$$

$$\hat{P}(x,s)=\frac{1}{v\rho }\sum _{k=0}^{\infty }\frac{(z-1{)}^k}{k!}$$

$$\displaystyle\times\frac{d^{k}}{dz^{k}}\left\{H_{0,1}^{1,0}\left[\rho^{2}z\left|\begin{array}[]{c l}-\\ (1,2)\end{array}\right.\right]\right\}_{z=1}$$

$$\displaystyle=\frac{1}{|x|}\sum_{k=0}^{\infty}\frac{1}{k!}\tau^{-k}s^{-k-1}H_{1,2}^{1,1}\left[\frac{x^{2}}{v^{2}}s^{2}\left|\begin{array}[]{c l}(0,1)\\ (1,2),(k,1)\end{array}\right.\right].$$

$$P(x,t)$$

$$P(x,t)\simeq \frac{1}{\sqrt{4\pi v^2\tau t}}\exp \left(-\frac{x^2}{4v^2\tau t}\right),$$

$$P(x,t)$$

$$=$$

$$⟨x^2(t)⟩={\mathcal{L}}^{-1}{\left[-\frac{{\partial }^2}{\partial k^2}\tilde{\hat{P}}(k,s)\right]|}_{k=0},$$

$$⟨x^2(t)⟩$$

$$E_{\alpha ,\beta }(z)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (\alpha k+\beta )},$$

$$⟨x^2(t)⟩\simeq \mathcal{D}\tau {\left(\frac{t}{\tau }\right)}^2,$$

$$⟨x^2(t)⟩\simeq 2\mathcal{D}t$$

$$⟨x^2(t)⟩$$

$$\overline{\stackrel{~}{\hat{P}}}(k_x,k_y,s)=\frac{1+s\tau -{\mathcal{D}}_x{k}_x^2\tilde{\hat{P}}(k_x,y=0,s)}{s(1+s\tau )+{\mathcal{D}}_y{k}_y^2},$$

$$\tilde{\hat{P}}(k_x,y,s)=\frac{1}{2\sqrt{{\mathcal{D}}_y}}\frac{\left[1+s\tau -{\mathcal{D}}_x{k}_x^2\tilde{\hat{P}}(k_x,y=0,s)\right]}{\sqrt{s(1+s\tau )}}$$

$$\times \exp \left(-\frac{\sqrt{s(1+s\tau )}}{\sqrt{{\mathcal{D}}_y}}|y|\right).$$

$$\tilde{\hat{P}}(k_x,y=0,s)=\frac{1}{2\sqrt{{\mathcal{D}}_y}}\frac{1+s\tau }{\sqrt{s(1+s\tau )}+\frac{{\mathcal{D}}_x}{2\sqrt{{\mathcal{D}}_y}}{k}_x^2},$$

$$\overline{\stackrel{~}{\hat{P}}}(k_x,k_y,s)$$

$$p_1(x,t)={\int }_{-\infty }^{\infty }P(x,y,t)dy.$$

$${\tilde{\hat{p}}}_1(k_x,s)=\overline{\stackrel{~}{\hat{P}}}(k_x,k_y=0,s)=\frac{1}{s}\frac{\sqrt{s(1+s\tau )}}{\sqrt{s(1+s\tau )}+\frac{{\mathcal{D}}_x}{2\sqrt{{\mathcal{D}}_y}}{k}_x^2},$$

Full text

1 Introduction

It is well known that any compact initial condition, evolving due to a parabolic (diffusion) equation, “spreads” instantly to infinity. This situation relates for example to the Fickian law of diffusion, or Fourier’s law of heat conduction, according equation $J=-\mathcal{D}\nabla P$, where $\mathcal{D}$ is the diffusion coefficient, $P=P(x,t)$ is the probability distribution function (PDF), while $J=J(x,t)$ is the probability flux. To overcome this unrealistic singularity with the infinite velocity of propagation, the so called telegrapher’s or Cattaneo equation has been introduced [1, 2]

[TABLE]

where $\tau$ is a characteristic time constant, $\mathcal{D}$ is the diffusion coefficient, which relates to a finite propagation velocity $v=\sqrt{\frac{\mathcal{D}}{\tau}}$. This equation is considered for various realizations of initial and boundary conditions. In the diffusion limit ($\tau\rightarrow 0$) of infinite-velocity propagation one recovers the standard diffusion equation. In the opposite limit $\tau\rightarrow\infty$, it becomes the wave equation with $\mathcal{D}/\tau=v^{2}$ being finite squared speed of the wave. The telegrapher’s equation was proposed by Kelvin and of Heaviside in electrodynamics theory and it was essentially employed in the heat transfer theory [3] and persistent random walk [1, 2, 3]. The standard Cattaneo equation has been solved analytically [4, 5], and further generalization on fractional order equations has been done by several authors [6, 7, 9, 11, 12, 13, 14, 15, 10, 8, 16].

Our main concern of the Cattaneo Eq. (1) in the comb geometry is motivated by ionic transport inside neuron dendrite structure. Recent experiments, together with numerical simulations, have investigated the calcium transport inside spiny dendrites to understand the role of calcium in signal transmission and neural plasticity. This issue is well reviewed in Refs. [17, 18]. Based on these experimental finding, different theoretical approaches have been developed to explore the transport properties of spiny dendrites. It is an active field of study [19, 20, 21, 22], and new experimental findings on calcium transport and reaction transport in neuroscience [22, 23] pose new questions to understand the impact of the geometry on calcium transport and reactions in spiny dendrites and the extension of various reaction-transport models to the case of subdiffusion. Recent experiments established a relation between the geometry of the dendrite spines and the subdiffusion observed in Refs. [24, 25]. It supports the application of the comb model, which is a convenient tool to explore analytically the anomalous transport in spiny dendrites in the framework of the continuous time random walk (CTRW) approach. Comb-like models can mimic ramified structures such as spiny dendrites of neuron cells [19, 20, 21] and can be used to describe the movement and binding dynamics of particles inside the dendritic spines. The corresponding process along the backbone is subdiffusiove with the power law evolution of the mean square displacement

[TABLE]

In this paper we consider finite-velocity diffusion on a comb, described by a comb model

[TABLE]

which represents the standard Cattaneo equation on a comb structure. The $\delta$-function means that the diffusion along the $x$-direction is allowed only at $y=0$ (the backbone). The particle moves in the backbone and can eventually be trapped in the fingers. The initial condition is given by

[TABLE]

and the boundary conditions for the PDF $P(x,y,t)$ and for $\frac{\partial}{\partial q}P(x,y,t)$, $q=\{x,y\}$ are set to zero at infinities, $x=\pm\infty$, $y=\pm\infty$. The diffusion coefficient along the $x$-direction is $\mathcal{D}_{x}\delta(y)$, with physical dimension $[\mathcal{D}_{x}\delta(y)]=\mathrm{m}^{2}/\mathrm{s}$, i.e., $[\mathcal{D}_{x}]=\mathrm{m}^{3}/\mathrm{s}$ ($[\delta(y)]=\mathrm{m}^{-1}$). The diffusion coefficient along the fingers is $\mathcal{D}_{y}$, $[\mathcal{D}_{y}]=\mathrm{m}^{2}/\mathrm{s}$. Correspondingly, the finite propagation velocities are $v_{x}=\mathcal{D}_{x}/[\tau\mathcal{D}_{y}]=\mathcal{D}_{x}/[v_{y}^{2}\tau^{2}]$ and $v_{y}=\sqrt{\mathcal{D}_{y}/\tau}$, where the $x$ component of the velocity accounts also the comb geometry, when transport along the backbone depends also on transport in fingers. Note also that the relation between the real three dimensional Laplace operator and the Laplace operator in the comb model (2) was established in Ref. [26].

2 Cattaneo equation

Prior concerning with the comb model (2) we discuss the properties of the telegrapher’s equation (1) with the initial conditions of Eq. (3) and natural (zero) boundary conditions at infinity. This equation was considered in Ref. [4], where an exact solution was obtained in the framework of Bessel functions. Here we give the solution in terms of the Fox $H$-function defined by inverse Mellin transform [27]:

[TABLE]

where

[TABLE]

and conditions of the integration are fulfilled111The conditions are $0\leq n\leq p$, $1\leq m\leq q$, $a_{i},b_{j}\in\mathrm{C}$, $A_{i},B_{j}\in\mathrm{R}^{+}$, $i=1,...,p$, $j=1,...,q$. Contour integration $\Omega$ starts at $c-\imath\infty$ and finishes at $c+\imath\infty$ separating the poles of the function $\Gamma(b_{j}+B_{j}s)$, $j=1,...,m$ with those of the function $\Gamma(1-a_{i}-A_{i}s)$, $i=1,...,n$..

After Fourier-Laplace transforms, one obtains the solution of Eq. (1) in the $(k,s)$-space

[TABLE]

Its inverse Fourier transform yields

[TABLE]

Using the variable change $z=1+(s\tau)^{-1}$ and $\rho=\frac{|x|}{v}s$, we present Eq. (10) in the compact form

[TABLE]

Now we apply Mellin transform to Eq. (11). Using a trick, we have

[TABLE]

This eventually yields the PDF in the Laplace space in the form of the Fox $H$-function

[TABLE]

Here we used the definition of the Fox $H$-function in Eq. (6). Expanding function $\hat{P}(x,s)$ in Taylor series about $z=1$, we have

[TABLE]

Application of the inverse Laplace transformation in Eq. (2) yields

[TABLE]

The asymptotic behavior of the solution in the case of $t/\tau\rightarrow\infty$ (i.e., $s\tau\rightarrow 0$) becomes Gaussian

[TABLE]

while the opposite case of $t/\tau\rightarrow 0$ (i.e., $s\tau\rightarrow\infty$), from Eq. (25), yields

[TABLE]

where we use the result given in [28] and properties of the Dirac delta function.

The MSD can be expressed in the form of the two parameter Mittag-Leffler function, as well. The MSD reads

[TABLE]

from where we find

[TABLE]

Here $E_{\alpha,\beta}(z)$ is the two parameter Mittag-Leffler function, defined by [30]

[TABLE]

where $z,\beta\in\mathrm{C}$, and $\Re(\alpha)>0$. Note that the two parameter Mittag-Leffler function (32) is a generalization of the one parameter Mittag-Leffler function $E_{\alpha}(z)=E_{\alpha,1}(z)$ and the exponential function $E_{1,1}(z)=e^{z}$.

For the short time limit $t/\tau\rightarrow 0$ the MSD corresponds to ballistic motion

[TABLE]

and then it changes to normal diffusion

[TABLE]

in the long time limit $t/\tau\rightarrow\infty$. These asymptotic diffusion and wave limits coincide exactly with results obtained in Ref. [4]. Here we note that the same result we obtain if one applies the exact solution (25), that is

[TABLE]

where we use the Mellin transform of the Fox $H$-function, in definition (8).

3 Solution of Cattaneo equation for a comb structure

Let us now analyze the comb Cattaneo equation (2) with the initial conditions (3) and zero boundary conditions at infinity. Fourier-Laplace transforming we find

[TABLE]

which yields after inverse Fourier transform in respect to $k_{y}$

[TABLE]

At this first step of finding a closed form of the solution, we need to find $\hat{P}(x,y=0,s)$. From Eq. (3) we have

[TABLE]

which eventually yields the closed form of the PDF in the Fourier-Laplace space

[TABLE]

Integrating Eq. (2) in respect to $y$ we obtain the marginal PDF

[TABLE]

Its Fourier-Laplace image results from Eq. (37) as follows222It is also known as the Montroll-Weiss equation [29].

[TABLE]

and correspondingly the Laplace image is

[TABLE]

The inverse Laplace transform can be found in the same way as it was done for the classical Cattaneo equation, cf. Eqs. (11)-(2). Rewriting Eq. (3) in the compact form by means of the variable change $z=1+(s\tau)^{-1}$ and $\rho=\sqrt{\frac{2\sqrt{\mathcal{D}_{y}\tau}}{\mathcal{D}_{x}}}|x|s^{1/2}$, we obtain $\hat{p}_{1}(x,s)$ in the form of the Fox H-function

[TABLE]

Its Taylor expansion about $z=1$ gives

[TABLE]

Now inverse Laplace transform yields the solution for the marginal PDF

[TABLE]

The asymptotic behavior of the PDF for $t/\tau\rightarrow\infty$ ($s\tau\rightarrow 0$) reads

[TABLE]

which is exactly the PDF for a diffusion equation on a comb [31]. Correspondingly, the short time solution for $t/\tau\rightarrow 0$ ($s\tau\rightarrow\infty$) becomes

[TABLE]

Comparing these limiting cases with corresponding limiting results in Eqs. (13) and (14), we conclude that the comb geometry affects strongly this wave-diffusion process. The wave dynamics is attenuated and fractional dynamics becomes dominant. It is also in a good qualitative agreement with experimental data on neocortical pyramidal neurons adapting with a time scale333We do not discuss this issue in the paper since it deserves a separate study., which is consistent with fractional order differentiation, such that the neuron’s firing rate is a fractional derivative of slowly varying stimulus parameters [32].

The fractional Cattaneo equation can be obtained as follows. Let us rewrite the Montroll-Weiss equation (38) as follows

[TABLE]

from where by inverse Fourier-Laplace transform we find the following generalized diffsuion equation

[TABLE]

where

[TABLE]

while $E_{\alpha,\beta}^{\delta}(z)$ is the three parameter Mittag-Leffler function [33]

[TABLE]

Here $\beta,\delta,z\in\mathrm{C}$, $\Re(\alpha)>0$, $(\delta)_{k}$ is the Pochhammer symbol $(\delta)_{0}=1$, $(\delta)_{k}=\frac{\Gamma(\delta+k)}{\Gamma(\delta)}$. Laplace transform of the Mittag-Leffler function reads

[TABLE]

Note also that $E_{\alpha,\beta}^{1}(z)=E_{\alpha,\beta}(z)$, and it is a special case of the Fox $H$-function [27]

[TABLE]

Eq. (57) can be rewritten as follows

[TABLE]

where

[TABLE]

is the regularized Prabhakar fractional derivative [34], and its Laplace transform reads

[TABLE]

This also results in the normal diffusion equation in the short time dynamics $t/\tau\rightarrow 0$,

[TABLE]

and it corresponds to the time fractional diffusion equation for $t/\tau\rightarrow\infty$,

[TABLE]

where ${{}_{C}}\mathcal{D}_{0+}^{1/2}$ is the Caputo fractional derivative [30]

[TABLE]

The MSD can be also estimated rigorously. From Eq. (56), we have

[TABLE]

Taking into account the series representation of the three parameter Mittag-Leffler function (59) we find that for the short time limit $t/\tau\rightarrow 0$ the MSD in Eq. (68) results in normal diffusion

[TABLE]

and in the long time limit $t/\tau\rightarrow\infty$ it changes to subdiffusion

[TABLE]

where we apply the asymptotic formula for the three parameter Mittag-Leffler function [35]

[TABLE]

with $z>1$, and $0<\alpha<2$.

It is worth noting that the MSD corresponds to the Mellin transform of the exact solution (50), which relates to Eq. (8) and yields

[TABLE]

4 Finite domain solution

In reality, the dendrite structure has a finite length along the backbone. In this case, the finite boundary conditions affect strongly the transient diffusion-subdiffusion process described by Eq. (62). One should also bear in mind that the spine-finger structure is finite as well. However, we suppose reasonably that the contaminant transport in this ramified structure is essentially slower than in the backbone, and the boundary conditions for the fingers are at $y=\pm\infty$. Therefore, now it is described by the marginal PDF in the framework of Eq. (62) with initial condition $p_{1}(x,t=0)=\delta(x)$ in the range $-L<x<L$ with absorbing boundary conditions $p_{1}(x=\pm L,t)=0$. It means that once a transporting particle reaches a boundary, it is instantly removed from the boundary.

We use the method of separation of variables $p_{1}(x,t)=X(x)T(t)$. Therefore, we find

[TABLE]

where $\lambda$ is a separation constant. From here we have a system of two equations

[TABLE]

Equation (73) is the eigenvalue problem with the boundary condition $X(x=\pm L)=0$, which yields a set of eigenfunctions $X_{n}$ with corresponding eigenvalues $\lambda_{n}$. Thus, the solution of Eq. (62) in the finite domain is

[TABLE]

where $T_{n}(t=0)=1\,,\forall n$. Eventually, after accounting for the initial condition $p_{1}(x,0)=\delta(x)$, the solution reads

[TABLE]

The solution of Eq. (72) can be found by the Laplace transform method. Thus, from relation (3), we have

[TABLE]

i.e.,

[TABLE]

The solution then becomes

[TABLE]

with $\lambda_{n}=\frac{\mathcal{D}_{x}}{2\sqrt{\mathcal{D}_{y}\tau}}\left[\frac{(2n+1)\pi}{2L}\right]^{2}$. We note that $T_{n}(t=0)=1$, since for $t=0$ only the first term with $j=0$ in the sum (77) survives. Thus, the finite domain solution reads

[TABLE]

Accounting relaxation in the finite boundaries, we also find the survival probability $S(t)=\int_{-L}^{L}p_{1}(x,t)\,dx$ which reads

[TABLE]

from where the first passage time PDF is

[TABLE]

The long time limit can be obtained by asymptotic expansion of the Mittag-Leffler function (69), which results in the following chain of transformations

[TABLE]

where we use the relation $E_{\alpha,\beta}(z)=zE_{\alpha,\alpha+\beta}(z)+\frac{1}{\Gamma(\beta)}$ [36]. From here, we find a power-law decay of the form $f(t)\simeq t^{-3/2}$, i.e.,

[TABLE]

where $\zeta(s,a)=\sum_{k=0}^{\infty}\frac{1}{(k+a)^{s}}$ is the Hurwitz zeta function [27].

The obtained result has a correct limit in the infinite domain $L\rightarrow\infty$. The solution (78) in the Laplace space is given by

[TABLE]

Taking the limit $L\rightarrow\infty$, i.e., $1/L\rightarrow 0$, the summation leads to integration which corresponds to the inverse Fourier transform with $k_{x}=\frac{(2n+1)\pi}{2L}$, from where we obtain an equivalent equation to Eq. (38) for the reduced PDF for the infinite domain case.

5 Summary

We have concerned with a finite velocity of a particle spreading in the framework of a Cattaneo equation in the comb geometry. This issue of anomalous diffusion with finite velocity in comb geometry was already established in numerical studies [13, 15, 37]. We presented exact solutions in the form of the Fox $H$-function series for both initial and boundary value problems. It describes kinetics, which is essentially complicated in comparison with anomalous diffusion described by comb fractional Fokker-Planck equation, obtained previously for a variety of realizations [38]. Fractional Cattaneo equation (62) describes a transient process of diffusion-subdiffusion with long time asymptotics of the MSD $\sim t^{1/2}$. For the dendrite transport described by Eq. (62), the boundary conditions are important. Various boundary value problems for the fractional Fokker-Planck equation have been already discussed in Ref. [39]. Here we considered a symmetrical boundaries with absorption, which is also important for the relaxation in many applications, including spiny dendrites and also for example in application of multichannel diffusion of hydrogen in solids [40].

Acknowledgements.

T.S. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG), project ME 1535/6-1 ”Random search processes, Lévy flights, and random walks on complex networks”. A.I. acknowledges the support from the Israel Science Foundation (ISF-931/16).