A unifying approach to first-passage time distributions in diffusing diffusivity and switching diffusion models

TL;DR

This paper introduces a unified theoretical framework to analyze first-passage time distributions in stochastic processes with evolving diffusivity, covering both continuous and discrete diffusivity models, with explicit formulas and solutions.

Contribution

It provides a comprehensive approach unifying diffusing diffusivity and switching diffusion models for first-passage time analysis, including general formulas and explicit solutions.

Findings

Derived moment-generating functions for integrated diffusivity

Provided explicit solutions for specific cases

Quantified diffusivity dynamics impact on first-passage times

Abstract

We propose a unifying theoretical framework for the analysis of first-passage time distributions in two important classes of stochastic processes in which the diffusivity of a particle evolves randomly in time. In the first class of "diffusing diffusivity" models, the diffusivity changes continuously via a prescribed stochastic equation. In turn, the diffusivity switches randomly between discrete values in the second class of "switching diffusion" models. For both cases, we quantify the impact of the diffusivity dynamics onto the first-passage time distribution of a particle via the moment-generating function of the integrated diffusivity. We provide general formulas and some explicit solutions for some particular cases of practical interest.

Equations

$$P(\mathbf{x},t|{\mathbf{x}}_0,D_0)=\sum _{n=1}^{\infty }{u}_n^{\ast }({\mathbf{x}}_0)u_n(\mathbf{x}){\rm Y}(t;{\lambda }_n|D_0),$$

$$T_t=\underset{0}{\overset{t}{\int }}d{t}^{\prime }{D}_{t^{\prime }},$$

$$P(\mathbf{x},t|{\mathbf{x}}_0)=\sum _{n=1}^{\infty }{u}_n^{\ast }({\mathbf{x}}_0)u_n(\mathbf{x}){\rm Y}(t;{\lambda }_n),$$

$$\rho (t|{\mathbf{x}}_0)=-\sum _{n=1}^{\infty }{u}_n^{\ast }({\mathbf{x}}_0){\partial }_t{\rm Y}(t;{\lambda }_n)\underset{\Omega }{\int }d\mathbf{x}{u}_n(\mathbf{x}).$$

$$d{D}_t=\mu (D_t,t)dt+\sigma (D_t,t)d{W}_t,$$

$$d{\mathbf{X}}_t=\sqrt{2D_t}d{\mathbf{W}}_t,$$

$${\rm Y}(t;\lambda |D_0,t_0)$$

$$\left({\partial }_{t_0}+\mu (D_0,t_0){\partial }_{D_0}+\frac{1}{2}{\sigma }^2(D_0,t_0){\partial }_{D_0}^2-\lambda D_0\right){\rm Y}(t;\lambda |D_0,t_0)=0,$$

$${\rm Y}(t;\lambda |D_0,t_0=t)=1,$$

$$\left.\biggl{(}\partial_{D_{0}}\Upsilon(t;\lambda|D_{0},t_{0})\biggr{)}\right|_{D_{0}=0}=0.$$

$$\mu (D,t)=\mu (D),\sigma (D,t)=\sigma (D).$$

$$\left({\partial }_t-\mu (D_0){\partial }_{D_0}-\frac{1}{2}{\sigma }^2(D_0){\partial }_{D_0}^2+\lambda D_0\right){\rm Y}(t;\lambda |D_0)=0,$$

$${\rm Y}(t;\lambda |D_0)=\sum _{n=0}^{\infty }\frac{(-1{)}^n}{n!}{\lambda }^n⟨T_t^n|D_0⟩$$

$$\left({\partial }_t-\mu (D_0){\partial }_{D_0}-\frac{1}{2}{\sigma }^2(D_0){\partial }_{D_0}^2\right)⟨T_t^n|D_0⟩=D_0⟨T_t^{n-1}|D_0⟩(n=1,2,\dots ),$$

$$\left({\partial }_t-\mu (D_0){\partial }_{D_0}-\frac{1}{2}{\sigma }^2(D_0){\partial }_{D_0}^2\right)⟨T_t|D_0⟩=D_0.$$

$$J_{\rm eq}(D)=\mu(D)p_{\rm eq}(D)-\partial_{D}\biggl{(}\frac{\sigma^{2}(D)}{2}\,p_{\rm eq}(D)\biggr{)}=0,$$

$${\rm Y}(t;\lambda )=\underset{0}{\overset{\infty }{\int }}d{D}_0{\rm Y}(t;\lambda |D_0)p_{\mathrm{eq}}(D_0).$$

$$\mu (D,t)=(\bar{D}-D)/\tau ,\sigma (D,t)=\sigma \sqrt{2D},$$

$${\rm Y}(t;\lambda |D_0)={\left(\frac{2\omega e^{-(\omega -1)t/(2\tau )}}{\omega +1+(\omega -1)e^{-\omega t/\tau }}\right)}^{\nu }$$

$$\times \exp \left(\frac{D_0(\omega +1)}{2{\sigma }^2\tau }\left(1-\frac{2\omega }{\omega +1+(\omega -1)e^{-\omega t/\tau }}\right)\right)$$

$${\rm Y}(t;\lambda )={\left(\frac{4\omega e^{-(\omega -1)t/(2\tau )}}{(\omega +1{)}^2-(\omega -1{)}^2{e}^{-\omega t/\tau }}\right)}^{\nu },$$

$$p_{\mathrm{eq}}(D)=\frac{{\nu }^{\nu }{D}^{\nu -1}}{\Gamma (\nu ){\bar{D}}^{\nu }}\exp (-\nu D/\bar{D}).$$

$$\mu (D,t)=0,\sigma (D,t)=\sigma ,$$

$$\biggl{(}\partial_{t}-\frac{\sigma^{2}}{2}\,\partial_{D_{0}}^{2}+\lambda D_{0}\biggr{)}\Upsilon(t;\lambda|D_{0})=0,$$

$$\left.\biggl{(}\partial_{D_{0}}\Upsilon(t;\lambda|D_{0})\biggr{)}\right|_{D_{0}=0}=\left.\biggl{(}\partial_{D_{0}}\Upsilon(t;\lambda|D_{0})\biggr{)}\right|_{D_{0}=D_{m}}=0.$$

$$v(D_{0})=C\biggl{[}\mathrm{Bi}^{\prime}(-\beta\gamma)\,\mathrm{Ai}\bigl{(}\alpha D_{0}-\beta\gamma\bigr{)}-\mathrm{Ai}^{\prime}(-\beta\gamma)\,\mathrm{Bi}\bigl{(}\alpha D_{0}-\beta\gamma\bigr{)}\biggr{]},$$

$$\underset{0}{\overset{D_m}{\int }}dD{v}^2(D)=1$$

$$\underset{0}{\overset{D_m}{\int }}dD{v}^2(D)=\frac{(\alpha D_m-\beta \gamma )v^2(D_m)+\beta \gamma v^2(0)}{\alpha },$$

$$C=\pi \sqrt{\alpha }{\left(\beta \gamma +(\alpha D_m-\beta \gamma ){\left(\frac{{\mathrm{Bi}}^{\prime }(-\beta \gamma )}{{\mathrm{Bi}}^{\prime }(\alpha D_m-\beta \gamma )}\right)}^2\right)}^{-1/2},$$

$$\mathrm{Ai}^{\prime}\bigl{(}-\beta\gamma\bigr{)}\mathrm{Bi}^{\prime}\bigl{(}\alpha D_{m}-\beta\gamma\bigr{)}-\mathrm{Ai}^{\prime}\bigl{(}\alpha D_{m}-\beta\gamma\bigr{)}\mathrm{Bi}^{\prime}\bigl{(}-\beta\gamma\bigr{)}=0.$$

$${\rm Y}(t;\lambda |D_0)=\sum _{k=0}^{\infty }{e}^{-{\sigma }^{2}t{\gamma }_k/2}{v}_k(D_0)\underset{0}{\overset{D_m}{\int }}dD{v}_k(D).$$

Full text

††: J. Phys. A: Math. Gen.

Keywords: diffusing diffusivity; switching diffusion; escape problem; first-passage time; diffusion-limited reaction

1 Introduction

An accurate description of biochemical reactions occuring in a heterogeneous, dynamically re-arranging intracellular environment is a long-standing problem [1, 2, 3, 4, 5, 6]. Various theoretical aspects of the underlying intracellular transport have been intensively studied over the past twenty years [7, 8]. In particular, different theoretical models have been proposed to account for molecular caging in the overcrowded cytoplasm [9, 10, 11, 12], viscoelastic properties of the cytoskeleton polymer network [13, 14, 15, 16, 17, 18], structural organization inside the cell [19], random diffusivity [20], and intermittent character of the motion [21, 22]. All these mechanisms affect the dynamics of molecules inside the cell, determine the statistics of their first-passage times (FPT) to the binding sites, and thus control the associated biochemical reactions.

Recently, we proposed a theoretical framework for investigating diffusion-limited reactions in dynamic heterogeneous media [23]. Modeling the effect of a rapidly re-arranging medium as random changes of the amplitude of thermal fluctuations felt locally by the tracer is based on the concept of diffusing diffusivity introduced by Chubynsky and Slater [24] and later explored by several authors [25, 26, 27, 28, 29]. Here, the diffusivity $D_{t}$ of the tracer is considered as a stochastic process, which is independent of the tracer’s position. For this annealed model, we derived a general expansion for the propagator $P(\bm{x},t|\bm{x}_{0},D_{0})$ of the tracer, i.e., the probability density of finding the tracer, started from $\bm{x}_{0}$ at time [math] with the initial diffusivity $D_{0}$, in a vicinity of a point $\bm{x}$ at time $t$:

[TABLE]

where the sum runs over all eigenvalues $\lambda_{n}$ and $L_{2}$-normalized eigenfunctions $u_{n}(\bm{x})$ of the Laplace operator $\Delta$ in a confined bounded medium $\Omega\subset{\mathbb{R}}^{d}$, $\Delta u_{n}+\lambda_{n}u_{n}=0$, with mixed Dirichlet/Neumann boundary conditions on the boundary $\partial\Omega$ of $\Omega$ [23]. As usual, the Dirichlet condition, $u_{n}(\bm{x})=0$ at $\bm{x}\in\Gamma\subset\partial\Omega$, accounts for a perfectly reactive sink $\Gamma$ on the boundary, whereas the Neumann condition $\partial_{\bm{n}_{\bm{x}}}u_{n}(\bm{x})=0$ at $\bm{x}\in\partial\Omega\backslash\Gamma$ describes an inert reflecting wall (an obstacle) on the remaining part of the boundary (here $\partial_{\bm{n}_{\bm{x}}}$ is the normal derivative at the boundary point $\bm{x}$ directed outward the domain).

While the structural and reactive properties of the medium are captured via the Laplacian eigenfunctions and eigenvalues [30], the function $\Upsilon(t;\lambda|D_{0})$ introduces the annealed disorder and couples it to the dynamics of the tracer. The function $\Upsilon(t;\lambda|D_{0})$ was shown to be the Laplace transform of the probability density function of the integrated diffusivity

[TABLE]

with the initial value $D_{0}$. For homogeneous diffusion with a constant diffusivity $D_{0}$, one gets $\Upsilon(t;\lambda)=\exp(-D_{0}t\lambda)$ and retrieves the standard spectral expansion of the propagator [31, 32]. If the initial diffusivity $D_{0}$ is chosen randomly from a prescribed distribution (e.g., the stationary distribution), the average of Eq. (1) yields the spectral expansion for the conventional propagator $P(\bm{x},t|\bm{x}_{0})$

[TABLE]

in which $\Upsilon(t;\lambda_{n})$ is the average of $\Upsilon(t;\lambda|D_{0})$, see below.

From the expansion (1), one easily gets the survival probability of the tracer, the macroscopic reaction rate and other quantities of interest. For instance, for a tracer started at $\bm{x}_{0}$, the probability density of the first-passage time to the binding site $\Gamma$ reads

[TABLE]

The impact of dynamic heterogeneities onto the first-passage time density is thus controlled by the function $\Upsilon(t;\lambda)$. When the diffusing diffusivity is modeled by a Feller process [33] (also known as the square root process or the Cox-Ingersoll-Ross process [34]), an explicit form of the function $\Upsilon(t;\lambda)$ was derived in [28] (see also Sec. 2.1). For this model, we obtained the asymptotic behavior of the probability density $\rho(t|\bm{x}_{0})$ and showed how the annealed disorder broadens the first-passage time distribution [23].

In the present paper, we further develop this theoretical approach by considering a general form of the stochastic equation for the diffusing diffusivity. Relying on the Feynman-Kac formula for the function $\Upsilon(t;\lambda|D_{0})$, we build a general framework for studying diffusion-limited reactions and the related first-passage time problems in the realm of diffusing diffusivity models (Sec. 2). In particular, we recall the main formulas for the Feller process and we derive new ones for the reflected Brownian motion on an interval. Moreover, we present similar results for another important class of models, in which a diffusing particle randomly switches between states with different diffusivities (Sec. 3). Such switching diffusion models are often employed to describe the dynamics in biological systems [22, 35, 36, 37, 38]. Quite naturally, switching diffusion models appear as discretized versions of diffusing diffusivity models. We formalize this connection by relating switching rates to the drift and volatility coefficients of the stochastic equation determining $D_{t}$. In this way, one gets a computationally efficient way to access the statistics of the first-passage time in both types of models.

2 Diffusing diffusivity models

We consider a particle diffusing in a $d$-dimensional dynamic heterogeneous environment, whose effect is modeled via a diffusing diffusivity $D_{t}$ which obeys a general stochastic equation in the Itô convention:

[TABLE]

subject to the initial condition $D_{t_{0}}=D_{0}$, where $W_{t}$ is the standard Wiener process, and functions $\mu(D,t)$ and $\sigma(D,t)$ represent drift and volatility of $D_{t}$, respectively. In turn, the position of the particle, $\bm{X}_{t}$, obeys another stochastic equation,

[TABLE]

in which $\bm{W}_{t}=(W_{t}^{1},\ldots,W_{t}^{d})$ is formed by independent Wiener processes $W_{t}^{i}$. In this basic setting that we employ throughout the paper, the particle undergoes locally isotropic diffusion driven by instantaneous interactions with the thermal bath (modeled by Gaussian noises $dW_{t}^{i}$) whose amplitude $\sqrt{2D_{t}}$ evolves with time. To account for inert impermeable walls or obstacles, a singular drift term should be added to the stochastic equation (6), see [39, 40] for technical details. More generally, one can include local anisotropy and external forces into Eq. (6) in a standard way [41]. These extensions would change the governing second-order differential operator and thus be incorporated via the modified eigenvalues and eigenfunctions, as for homogeneous diffusion. In contrast, the effect of rapid re-arrangements of the medium on larger length scales is captured by the diffusing diffusivity $D_{t}$, which is considered to be independent of local thermal noises $dW_{t}^{i}$. Most importantly, the stochastic equation (5) does not depend on the position $\bm{X}_{t}$ of the tracer. As a consequence, one can study separately the dynamics of the diffusivity and then subordinate the dynamics of the tracer according to Eqs. (1, 4).

For this purpose, one needs to evaluate the moment-generating function $\Upsilon(t;\lambda|D_{0},t_{0})$, which is the Laplace transform of the probability density $Q(t,T|D_{0},t_{0})$ of the integrated diffusivity $T_{t}$ defined in Eq. (2):

[TABLE]

given that the initial diffusivity at $t_{0}$ is $D_{0}$. This function satisfies the Feynman-Kac formula

[TABLE]

subject to the terminal condition

[TABLE]

an appropriate boundary condition at $D_{0}=0$, and a regularity condition $\Upsilon(t;\lambda|D_{0},t_{0})\to 0$ as $D_{0}\to\infty$. The boundary condition at $D_{0}=0$ should ensure that the diffusivity remains nonnegative. Following the discussion in [28, 23], we impose the no-flux boundary condition to maintain the normalization of the probability density for diffusivity. For the backward equation, this condition reads

[TABLE]

We note that the terminal condition (9) postulates the average over the diffusivity at time $t$. If one was interesting in knowing the value of diffusivity at time $t$, $D_{t}=D$, the terminal condition would be replaced by $\Upsilon(t;\lambda|D_{0},t_{0}=t)=\delta(D-D_{0})$. In particular, substituting such $\Upsilon(t;\lambda|D_{0},t_{0})$ into the spectral expansion Eq. (1) would yield the full propagator $P(\bm{x},D,t|\bm{x}_{0},D_{0},t_{0})$ characterizing both the position and the diffusivity of the particle. However, we do not consider this extension in the paper.

In the remaining part of the paper, we focus on the case of diffusing diffusivity that is homogeneous in time, i.e., the drift and volatility coefficients are time-independent:

[TABLE]

In this case, the solution of Eq. (8) depends on the difference $t-t_{0}$, i.e., $\Upsilon(t;\lambda|D_{0},t_{0})=\Upsilon(t-t_{0};\lambda|D_{0},0)$, so that one can replace $\partial_{t_{0}}$ by $-\partial_{t}$ and then set $t_{0}=0$:

[TABLE]

subject to the initial condition $\Upsilon(t=0;\lambda|D_{0})=1$ and the same boundary condition (10) and regularity condition (note that we omitted $t_{0}=0$ in $\Upsilon(t;\lambda|D_{0})$). This equation for the function $\Upsilon(t;\lambda|D_{0})$ can also be used to derive equations for the moments of the integrated diffusivity $T_{t}$. In a standard way, substituting the expansion

[TABLE]

into Eq. (12) and grouping the terms of the same order in $\lambda$ yield a set of equations

[TABLE]

subject to the initial condition $\langle T_{0}^{n}|D_{0}\rangle=0$ and former boundary conditions. In particular, the mean integrated diffusivity obeys

[TABLE]

When there exists a unique equilibrium distribution of diffusivity, $p_{\rm eq}(D)$, at which the probability flux of the associated forward Fokker-Planck equation (with $\lambda=0$) vanishes, i.e.

[TABLE]

it is convenient to average over the initial diffusivity $D_{0}$ drawn from this equilibrium density:

[TABLE]

2.1 Example: a Feller process

In [28, 23], we studied in detail the diffusing diffusivity modeled by a Feller process, for which

[TABLE]

with three parameters: the mean diffusivity $\bar{D}$, the relaxation time scale $\tau$, and the amplitude of diffusivity fluctuations $\sigma$. In particular, we obtained111 Three misprints were found in [28] in inline equations after Eq. (6), compare them with our corrected Eq. (19). These misprints did not affect the remaining content of Ref. [28].

[TABLE]

and

[TABLE]

with $\omega=\sqrt{1+4\sigma^{2}\tau^{2}\lambda}$, $\nu=\bar{D}/(\tau\sigma^{2})$, and the equilibrium diffusivity is known to follow the Gamma distribution:

[TABLE]

The first-passage and extreme value properties of the Feller process itself were studied earlier in [42, 43, 44].

2.2 Example: reflected Brownian motion

Chubynsky and Slater first introduced the diffusing diffusivity qualitatively as reflected Brownian motion on the positive half-line [24]. To avoid an unlimited growth of diffusivity, it is more natural to consider Brownian motion on an interval $(0,D_{m})$ with two reflecting endpoints. We explore this case with

[TABLE]

so that Eq. (12) becomes

[TABLE]

subject to the initial condition $\Upsilon(t=0;\lambda|D_{0})=1$ and two boundary conditions:

[TABLE]

The amplitude of fluctuations, $\sigma$, strongly affects the diffusivity dynamics. In the limit $\sigma\to 0$, fluctuations are suppressed, and one deals with a constant initial diffusivity $D_{0}$. In the opposite limit $\sigma\to\infty$, the diffusivity switches so rapidly between different values in $(0,D_{m})$ that its behavior resembles a constant mean diffusivity $D_{m}/2$.

The equation (23) admits a standard spectral solution in terms of eigenvalues and eigenfunctions of the associated differential operator $\mathcal{L}=\partial_{D_{0}}^{2}-\alpha^{3}D_{0}$ (with $\alpha^{3}=2\lambda/\sigma^{2}$). One can search for such an eigenpair $\gamma$ and $v(D_{0})$ in the form

[TABLE]

where $\mathrm{Ai}(z)$ and $\mathrm{Bi}(z)$ are two linearly independent Airy functions, prime denotes the derivative, $\beta=1/\alpha^{2}$, and $C$ is a normalization constant which is fixed by imposing the $L_{2}$-normalization of the eigenfunction:

[TABLE]

(see Refs. [45, 46] for a more detailed analysis of a similar problem). This integral can be evaluated using the Airy equation and the reflecting boundary conditions:

[TABLE]

from which $C$ can be expressed as

[TABLE]

where we used the Wronskian of Airy functions and reflected boundary conditions (note also that $v(0)=C/\pi$ and $v(D_{m})=C\mathrm{Bi}^{\prime}(-\beta\gamma)/(\pi\mathrm{Bi}^{\prime}(\alpha D_{m}-\beta\gamma))$).

The form (25) already satisfies the reflecting boundary condition at $D_{0}=0$. The eigenvalue $\gamma$ is determined from the second boundary condition at $D_{0}=D_{m}$ that implies

[TABLE]

As $-\alpha^{3}D_{0}$ is a bounded perturbation of the double derivative operator on an interval, the spectrum of the operator $\mathcal{L}$ is discrete, i.e., there are infinitely many solutions $\gamma_{k}$ of Eq. (29) that we enumerate by index $k=0,1,2,\ldots$. The associated eigenfunctions $v_{k}$ form a complete orthonormal basis in $L_{2}(0,D_{m})$. As a consequence, the solution of Eq. (23) can be decomposed on this basis as

[TABLE]

If the initial diffusivity $D_{0}$ is drawn from the equilibrium density $p_{\rm eq}$ (which is uniform in this setting), one gets

[TABLE]

From the moment-generating function, one can compute the moments of the integrated diffusivity. One can either perform the asymptotic analysis of Eq. (30) as $\lambda\to 0$, or solve directly Eq. (15) with $\mu=0$, subject to the reflecting boundary conditions at [math] and $D_{m}$. In the latter case, an expansion of $\langle T_{t}|D_{0}\rangle$ over the complete basis of cosine functions on $(0,D_{m})$ leads to

[TABLE]

In the short-time limit, one retrieves $\langle T_{t}|D_{0}\rangle\simeq D_{0}t$, whereas in the long-time limit, one gets

[TABLE]

with the expected dominant behavior $D_{m}t/2$. Higher-order moments obeying Eqs. (14) can be found in the same way.

In sharp contrast to the fully explicit solution (19) for the Feller process, the solution (30) requires a numerical computation of eigenvalues $\gamma_{k}$ which depend implicitly on the parameter $\lambda$. As the propagator $P(\bm{x},t|\bm{x}_{0})$ in Eq. (1) also involves a spectral decomposition over the Laplacian eigenvalues $\lambda_{n}$ in a bounded domain, the eigenvalues $\gamma_{k}$ should be evaluated for each $\lambda_{n}$ that makes this solution computationally demanding and impractical. At the same time, this solution allows one to analyze the asymptotic behavior of $\Upsilon(t;\lambda)$ and all related quantities as we briefly illustrate below.

For small $\lambda$, the term $-\alpha^{3}D_{0}$ can be considered as a small perturbation of the double derivative in the operator $\mathcal{L}=\partial^{2}_{D_{0}}-\alpha^{3}D_{0}$ so that the eigenvalues and eigenfunctions of $\mathcal{L}$ are close to that of the double derivative operator. The perturbation theory yields thus

[TABLE]

As expected, the correction terms are small for the modes with $k=1,2,\ldots$, whereas the correction term is dominant for the constant mode with $k=0$. In this limit, one gets

[TABLE]

given that the contribution of other terms is small because the (unperturbed) eigenfunctions $\cos(\pi kD/D_{m})$ are orthogonal to $1$. This analysis is also applicable in the limit $\sigma\to\infty$, in which fluctuations are so strong that the model is reduced to homogeneous diffusion with the mean diffusivity $D_{m}/2$.

In the opposite limit of large $\lambda$, one deals with large $\alpha$ and small $\beta$ in Eq. (29) so that the function \mathrm{Bi}^{\prime}\bigl{(}\alpha D_{m}-\beta\gamma\bigr{)} is exponentially large, whereas the function \mathrm{Ai}^{\prime}\bigl{(}\alpha D_{m}-\beta\gamma\bigr{)} is exponentially small. As a consequence, zeros of Eq. (29) are very close to the zeros of \mathrm{Ai}^{\prime}\bigl{(}-\beta\gamma\bigr{)}, i.e.,

[TABLE]

where $a^{\prime}_{k}$ are the zeros of the derivative of the Airy function (e.g., $a^{\prime}_{0}\approx-1.0188$). We conclude that both $\Upsilon(t;\lambda|D_{0})$ and $\Upsilon(t;\lambda)$ decay with $\lambda$ in a stretched-exponential way.

One can see that the ratio $\sigma^{2}/D_{m}^{3}$, setting the borderline between two asymptotic limits (34, 36), introduces a characteristic length of dynamic disorder, $\sqrt{D_{m}^{3}/\sigma^{2}}$. This length scale is compared in Eqs. (1, 4) to the diffusion length $\sqrt{D_{m}t}$ and to the geometric length scales of the reactive medium determined by the eigenvalues $\lambda_{n}^{-1/2}$. In particular, various asymptotic limits of the first-passage time density can be deduced, in analogy with the results presented in [23] for the case of the diffusivity modeled by a Feller process.

It is instructive to look at the limit $\sigma=0$, in which the diffusivity does not fluctuate, so that $\Upsilon(t;\lambda|D_{0})=e^{-D_{0}\lambda t}$, where $D_{0}$ is the initial diffusivity. If this diffusivity is randomly chosen from the equilibrium distribution, one gets

[TABLE]

As a consequence, when either $t$ or $\lambda$ goes to infinity, the function $\Upsilon(t;\lambda)$ decays slowly (as a power law). This slow decay is a consequence of the superstatistical description: the average over $D_{0}$ includes the contribution from particles with arbitrarily small diffusivities. This is drastically different from the stretchted-exponential decay with respect to $\lambda$ and from the exponential decay with respect to $t$ in the presence of fluctuations: even though small diffusivities are still accessible, it is unlikely that the particle keeps such a small diffusivity for a long time.

We also note that the solution of a more general problem of reflected Brownian motion on a shifted interval $(D_{1},D_{1}+D_{m})$ (with $D_{1}>0$) can be easily reduced to our solution by shifting the diffusivity $D_{t}$, i.e., by considering $\hat{D}_{t}=D_{t}+D_{1}$, where $D_{t}$ is modeled by reflected Brownian motion on $(0,D_{m})$ as before. The shift by $D_{1}$ leads to a constant term $D_{1}t$ in the integrated diffusivity so that $\hat{\Upsilon}(t;\lambda)=e^{-D_{1}t\lambda}\Upsilon(t;\lambda)$. We note that the explicit factor $e^{-D_{1}t\lambda}$ provides the dominant contribution to the decrease of the function $\hat{\Upsilon}(t;\lambda)$ as compared to $\Upsilon(t;\lambda)$ in the limit $\lambda\to\infty$.

2.3 Moments of the position

In [28], the propagator $P(x,t|x_{0})$ for diffusion on the line (without boundary) was expressed in terms of the moment-generating function $\Upsilon(t;\lambda)$ for the Feller process:

[TABLE]

The subordination argument [23, 27] supports this relation for any model of diffusing diffusivity. This relation provides thus an additional interpretation of $\Upsilon(t;\lambda)$ as the characteristic function of the one-dimensional displacement on a line. In particular, one can easily evaluate the moments of the displacement, e.g.,

[TABLE]

As expected, odd moments vanish due to the symmetry of thermal noise $dW_{t}$ (and independently of the diffusivity model), whereas even moments can be expressed through the moments of the integrated diffusivity

[TABLE]

This is an extension of the basic relation for the moments of the homogeneous Gaussian diffusion, for which $\langle T_{t}^{n}\rangle=(D_{0}t)^{n}$. From this relation, one easily gets the kurtosis, as well as the non-Gaussian parameter, $\gamma(t)=\langle X_{t}^{4}\rangle/(3\langle X_{t}^{2}\rangle^{2})-1$. Note that the same relation holds for the moments with a prescribed initial diffusivity $D_{0}$ which can be found by solving Eqs. (14):

[TABLE]

In contrast, there is no such a direct relation between the moments $\langle T_{t}^{n}\rangle$ and $\langle X_{t}^{2n}\rangle$ for restricted diffusion.

3 Switching diffusion model

A switching diffusion model, in which a diffusing particle randomly switches between internal states with distinct diffusivities, can be considered as a discrete version of diffusing diffusivity models. According to the subordination argument [23], it is enough to obtain the propagator for one-dimensional switching diffusion on a line, see Eq. (38), whereas its Fourier transform yields the function $\Upsilon(t;\lambda)$ and thus accesses general FPT problems in arbitrary confined reactive domains in ${\mathbb{R}}^{d}$ via Eqs. (1, 4). We consider such a model with $J$ states which are characterized by a set of diffusion coefficients $D_{i}$ and switching rates $k_{ij}$ (a rigorous mathematical formulation of switching models and some their properties can be found in [22, 37, 47, 48]). We introduce the probability density $P_{i,i_{0}}(x,t|x_{0},t_{0})$ of finding the particle in a vicinity of $x$ in the state $i$ at time $t$, given that it was started from $x_{0}$ in the state $i_{0}$ at time $t_{0}$. This propagator satisfies the forward Fokker-Planck equation

[TABLE]

subject to the initial condition: $P_{i,i_{0}}(x,t=t_{0}|x_{0},t_{0})=\delta_{i,i_{0}}\delta(x-x_{0})$, where we defined $k_{ii}=-\sum\limits_{j\neq i}k_{ij}$. The first term on the right-hand side describes diffusion (with diffusivity $D_{i}$), while the second term accounts for switching between different states. In this class of switching models, the dynamics in each internal state is governed by the same differential operator and differs only by its diffusivity. This is the crucial property that will allow for getting the propagator $P(\bm{x},t|\bm{x}_{0})$ in Eq. (1) for this model. Such an extension is not directly applicable to other intermittent processes, in which the governing operator changes between states. Moreover, modifications are needed even in the case when the operator remains the same (e.g., the Laplace operator) but the boundary condition changes between states (see, e.g., a two-state model developed in [38], in which the particle reacts with the target only when it is in an “active” state). Similarly, the models of surface-mediated diffusion [49, 50, 51, 52, 53] are not considered here as their switching mechanisms are different.

The Fourier transform reduces the partial differential equations (45) to a set of linear ordinary differential equations that can be solved in a matrix form, from which

[TABLE]

where ${\bf D}$ is the diagonal $J\times J$ matrix of diffusivities, ${\bf D}_{ij}=\delta_{ij}D_{i}$, and ${\bf K}$ is the matrix of switching rates, ${\bf K}_{ij}=k_{ij}$. If $p_{i}$ denotes the probability of starting in the initial state $i$, the marginal propagator averaged over the initial and arrival states reads

[TABLE]

where

[TABLE]

Note that the function $\Upsilon(t;\lambda|D_{0})$ admits the same form, with $p_{j}$ being equal to [math] for all states, except for the state with $D_{0}$ (for which $p_{j}=1$).

Using this relation, one can access the propagator $P(\bm{x},t|\bm{x}_{0})$ and the first-passage time density $\rho(t|x_{0})$ in a general confining domain according to Eqs. (1, 4). We recall that the matrix exponential function in Eq. (48) can be evaluated via diagonalization of the matrix $\lambda{\bf D}-{\bf K}^{\dagger}$. A fully explicit solution can be obtained for the two-state switching model (see, e.g., [28, 54]):

[TABLE]

where $\bar{D}=p_{1}D_{1}+p_{2}D_{2}$ and

[TABLE]

This rigorous result refines a former discussion of the two-state noise in [55]. For a larger number of states, formulas rapidly become too cumbersome and impractical. In contrast, a numerical computation of the function $\Upsilon(t;\lambda)$ via the matrix form (48) remains efficient even for relatively large number of states (up to few thousand).

3.1 Relation to diffusing diffusivity models

In this subsection, we discuss how discretized versions of diffusing diffusivity models with time-independent coefficients $\mu(D)$ and $\sigma(D)$ are related to switching diffusion models. In fact, continuously varying diffusivity $D_{t}$ can be replaced by a set of discrete values $D_{i}=i\,{\epsilon}$ ($i=1,2,\ldots$), with a discretization step ${\epsilon}$. Variations of $D_{t}$ can thus be seen as switching between neighboring states. We briefly discuss two equivalent approaches to formalize this connection.

In the first approach, the switching rates are determined by discretizing the forward Fokker-Planck equation for the probability density $p(D,t|D_{0},t_{0})$ of the diffusivity $D_{t}$ (in the Itô convention)

[TABLE]

subject to the initial condition $p(D,t=t_{0}|D_{0},t_{0})=\delta(D-D_{0})$. The discretization of the right-hand side of this equation with a diffusivity step ${\epsilon}$ reads

[TABLE]

In other words, the original PDE is approximated by a set of linear ordinary differential equations. These discretized equations for $p(i{\epsilon},t|D_{0},t_{0})$ can be seen as a switching model with multiple states $D_{i}=i{\epsilon}$ ($i=1,2,\ldots$) and the switching rates

[TABLE]

and zero otherwise, where we used the shortcut notations $\sigma_{i}=\sigma(i{\epsilon})$ and $\mu_{i}=\mu(i{\epsilon})$. As the “first” equation for $p({\epsilon},t)$ involves the term $p(0,t)$, one needs to close this system by accounting for the reflecting boundary condition at $D=0$:

[TABLE]

Expressing $p(0,t)$ in terms of $p({\epsilon},t)$ leads to a slight modification of the first diagonal element: $k_{1,1}=-\sigma^{2}_{1}/(2{\epsilon}^{2})-\mu_{1}/{\epsilon}$.

As a numerical solution of an infinitely-dimensional system of equations is not feasible, one needs to truncate the original problem by imposing an additional reflecting boundary condition at some truncation level $D_{m}$. As in the case of $D=0$, this boundary condition changes the coefficient $k_{J,J}$ in the “last” equation for $p(D_{m},t)$, with $D_{m}=J{\epsilon}$. However, a simpler and more consistent way of closing the system of equations is to require $\sum\nolimits_{j}k_{ij}=0$ for $i=J$. This is a detailed balance condition for the matrix of switching rates, which is already satisfied for all rates from Eq. (52) with $i=1,\ldots,J-1$. Combining these relations, one can write the set of discretized equations in a matrix form as

[TABLE]

with the $J\times J$ three-diagonal matrix

[TABLE]

whereas ${\bf D}_{ij}=\delta_{ij}i{\epsilon}$. In other words, we identified the matrices ${\bf D}$ and ${\bf K}$ determining a switching diffusion model that is a discrete approximation of the diffusing diffusivity model with coefficients $\mu(D)$ and $\sigma(D)$. The relation (48) expresses the function $\Upsilon(t;\lambda)$ for this switching model.

Alternatively, one could directly discretize the backward equation (8):

[TABLE]

(here for convenience we adopted another discretization scheme for the first derivative). The solution of this discretized equation with the terminal condition $\Upsilon(t;\lambda|D_{0},t_{0}=t)=1$ reads

[TABLE]

with the matrices ${\bf D}$ and ${\bf K}$ defined above. Averaging this solution over the initial states chosen with probabilities $p_{j}$ and setting $t_{0}=0$ yield

[TABLE]

which is just a transposed re-writing of Eq. (48).

3.2 Numerical illustrations

Figure 1 illustrates a comparison between a diffusing diffusivity model and its approximation by switching diffusion. We set the drift and the volatility terms according to Eq. (18) for the Feller process, with $\bar{D}=1$, $\sigma=1$, and $\tau=1$ (arbitrary units). On one hand, the function $\Upsilon(t;\lambda)$ is computed via explicit analytical solution (20). On the other hand, a discrete approximation of this process by switching diffusion allows one to compute the function $\Upsilon(t;\lambda)$ by Eq. (48). This computation depends on the discretization step ${\epsilon}$ and the truncation threshold $D_{m}$. We set $D_{m}=10$ and checked that further increase of this value does not almost affect the computation. This is not surprising given that the equilibrium distribution of diffusivities, Eq. (21), decays exponentially fast for the Feller process. The discretization step has also relatively weak impact on the solution, if it is small enough (see Fig. 1). We emphasize, however, the quality of the approximation depends in general on the chosen model and its parameters. For instance, setting $\tau=10$ (while keeping $\bar{D}=1$ and $\sigma=1$) yields $\nu=0.1$ in the Feller model and thus an integrable but divergent at $D=0$ equilibrium density in Eq. (21). As a consequence, a much finer discretization is needed to accurately capture the behavior of this density near zero and thus to get an accurate representation via a switching diffusion model. In general, one needs to undertake the convergence analysis or at least to compute $\Upsilon(t;\lambda)$ with various discretization steps ${\epsilon}$ to check its convergence.

Left panels of Fig. 2 show the behavior of the function $\Upsilon(t;\lambda)$ for a diffusing diffusivity modeled by reflected Brownian motion on $(0,D_{m})$. Although the exact solution is provided in Sec. 2.2, it is much faster and easier to use the approximate solution via the switching diffusion model. As expected, the function $\Upsilon(t;\lambda)$ approaches $1$ as $t\to 0$ or $\lambda\to 0$. In turn, when either of these variables is getting large, $\Upsilon(t;\lambda)$ decreases. According to Eq. (37), the decay is slow for the case without diffusivity dynamics ($\sigma=0$), see Fig. 2(a). In turn, much faster decay is observed for other cases with $\sigma>0$, in agreement with the asymptotic analysis of Sec. 2.2.

To get a closer look into the behavior of $\Upsilon(t;\lambda)$, it is convenient to plot $-\ln(\Upsilon(t;\lambda))/(t\lambda)$ as a function of $t$. At small $t$, one has $T_{t}\approx D_{0}t$ so that $\Upsilon(t;\lambda)\simeq\langle e^{-t\lambda D_{0}}\rangle\simeq 1-t\lambda\langle D_{0}\rangle$, where $\langle D_{0}\rangle$ is the mean initial diffusivity, which is equal to $D_{m}/2$ in this model. As a consequence, the ratio $-\ln(\Upsilon(t;\lambda))/(t\lambda)$ approaches $D_{m}/2$ as $t\to 0$. The opposite limit $t\to\infty$ is less universal: for instance, Eq. (37) exhibits $\ln(t)/t$ decay, whereas Eq. (31) leads to a constant. Right panels of Fig. 2 illustrate this behavior. All shown curves start from the mean diffusivity $D_{m}/2=0.5$ at $t=0$. As $t$ grows, all curves decrease but the speed of decrease depends on $\lambda$ and $\sigma$. For $\sigma=0$, the ratio vanishes as $t\to\infty$ whereas it reaches a nonzero limit $\sigma^{2}\gamma_{0}(\lambda)/(2\lambda)$ for $\sigma>0$, where $\gamma_{0}(\lambda)$ is the smallest eigenvalue of the operator $\partial^{2}_{D_{0}}-(2\lambda/\sigma^{2})D_{0}$, see Sec. 2.2. As a consequence, this limit changes from $D_{m}/2$ for small $\lambda$ to $|a^{\prime}_{0}|(\sigma^{2}/2)^{1/3}\lambda^{-1/3}$ for large $\lambda$. One can see that the range of variations of $-\ln(\Upsilon(t;\lambda))/(t\lambda)$ is getting narrower as $\sigma$ increases. Indeed, strong fluctuations rapidly mix all diffusivities in $(0,D_{m})$, restoring the behavior with the mean diffusivity $D_{m}/2$, as at short times.

4 Conclusion

We formulated a unifying approach for studying first-passage time distributions and related diffusion-limited reactions in the realm of diffusing diffusivity and switching diffusion models. In both cases, the dynamics of randomly changing diffusivity $D_{t}$ is assumed to be independent from the particle’s position, so that the subordination argument yields a general spectral expansion for the propagator and the first-passage time probability density. The key element coupling the stochastic diffusivity to the motion of the particle is the moment-generating function $\Upsilon(t;\lambda|D_{0})$ of the integrated diffusivity.

In diffusing diffusivity models, continuous changes of $D_{t}$ are governed by a stochastic differential equation, and the function $\Upsilon(t;\lambda|D_{0})$ can be calculated by using the Feynman-Kac formula. We illustrated this formalism for the case when the diffusivity is modeled by reflected Brownian motion on an interval with reflecting endpoints. In turn, when the diffusivity randomly switches between discrete values, we derived a matrix representation of the function $\Upsilon(t;\lambda|D_{0})$ involving the matrix of switching rates. We also formalized the connection between these two classes of models by relating the coefficients of the stochastic differential equation, $\mu(D)$ and $\sigma(D)$, to the switching rates. With the help of this formalism, one can extend former results on diffusive search problems [2, 5, 32, 56, 57, 58] to heterogeneous diffusion, compute the related first-passage time distributions [3, 59, 60, 61] and investigate diffusion-limited reactions in heterogeneous media.

While continuously changing diffusivity may represent the effect of rapidly re-arranging medium onto the motion of particles, discrete changes of the diffusivity can mimic switching between conformational states of a polymer or reversible binding of the diffusing molecule to other constituents (static or mobile) of the medium. In particular, the state with zero diffusivity can incorporate trapping events. While former studies involving stochastic diffusivity were focused on a specific choice of the Feller process (which includes as a particular case the square of the Ornstein-Uhlenbeck process used in [25, 26, 27]), the general formalism of the present paper opens the door to study a very broad class of various processes in a unified way. As the microscopic theory expressing the impact of rapidly re-arranging media onto the particle’s dynamics in terms of an appropriate diffusing diffusivity model is still missing, the possibility of dealing with a broad range of “candidate processes” is particularly valuable for future research.