Yann Lanoiselée
email: y.lanoiselee@bham.ac.uk

Denis S Grebenkov
email: denis.grebenkov@polytechnique.edu

Non-Gaussian diffusion of mixed origins

Keywords: non-Gaussian diffusion, heterogeneous media, diffusing diffusivity, anomalous diffusion, first-passage times

The properties of diffusion processes are drastically affected by heterogeneities of the medium that can induce non-Gaussian behavior of the propagator in contrast with the idealized realm of Brownian motion. In this paper we analyze the diffusion propagator when distinct origins of heterogeneity (e.g. time-fractional diffusion, diffusing diffusivity, distributed diffusivity across a population) are combined. These combinations allow one to describe new classes of strongly heterogeneous processes relevant to biology. Based on a combined subordination technique, we obtain the exact propagator for different instructive examples. This approach is then used to calculate analytically the first-passage time statistics (on half-real line and in any bounded domain) for a particle undergoing non-Gaussian diffusion of mixed origins.

Introduction

The description of diffusion processes in heterogeneous media is fundamental for characterizing physical and biological diffusion-controlled systems. In the last decade, many experiments in complex media have shown evidences that heterogeneities have a drastic effect onto the diffusive properties of particles. In particular, the distribution of displacements of a particle can significantly deviate from the conventional Gaussian distribution [1][2][3][4][START_REF] Bertrand | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Moschakis | [END_REF][7][8][START_REF] Wang | Proc. Nat. Acad. Sci[END_REF][START_REF] Wang | [END_REF][11][12][13]. The non-Gaussian character, observed at short times, can either disappear at long times or not. In the former case, the particle can be assumed to explore different diffusivities along its path that induce non-Gaussian displacements. In such theoretical models, when the diffusivity is fluctuating around its mean [14][15][16][17][18][19][20], the Gaussian regime is recovered when the trajectory duration is much longer than the meanreverting time of diffusivity [21,22]. In another group of models, known as annealed transit time model (ATTM), the particle experiences a medium where the diffusion coefficient does not evolve continuously but changes at random times [23].

When diffusion remains non-Gaussian even at long times, several explanations can be at play. When the distribution of periods with a constant diffusivity in the ATTM model has a heavy tail with infinite mean, the particle never explores the whole distribution of diffusivities. Therefore the dynamics exhibits weak ergodicity breaking [23,24]. Another model that displays both non-Gaussian distribution and weak ergodicity breaking is the time-fractional diffusion when the particle dynamics is so slow that it never explores the whole phase space. The non-Gaussian behavior can also originate from intertrajectory fluctuations when trajectories with different diffusion coefficients are used to build up an empirical probability density [26][27][28]. This situation can typically happen in cell biology when inter-cell fluctuations are unavoidable. Another occurrence of non-Gaussian distributions arises when a Gaussian process is superimposed with additional non-Gaussian noise. A priori, all these mechanisms are not mutually exclusive and can combine to produce very complicated types of motion.

In this article, we investigate the probability density function of particle displacements subject to non-Gaussian diffusion of mixed origins when several mechanisms contribute.

Previously, mixed diffusion processes has already been discussed for the particular cases of Continuous Time Random Walk on a fractal structure [29] and subordinated fractional Brownian motion [30,31]. We will take advantage of the fact that many non-Gaussian diffusion processes can be mathematically described via subordination [32][33][34][35][36][37]. Then we successively apply the subordination technique in order to describe the combined effect of different mechanisms causing heterogeneous non-Gaussian diffusion. The path properties of the particle remain the same as in the case of Brownian motion but the "speed" at which the particle travels along the path is disorder-driven. An equivalent way of interpreting the subordination is that the integrated diffusivity is not deterministic but randomly increasing.

Sec. 2, summarizes in a systematic way several formerly known physical mechanisms leading to non-Gaussian diffusion in terms of subordination integrals and describes the particular forms of the probability density function of their integrated diffusivity. In Sec.

3, we introduce a model of integrated diffusivity which can either converge to normal diffusion at long time or not. We also discuss its relation to other random diffusivity models. In Sec. 4, using the properties of subordination, we obtain the probability density function of particle displacements, subject to several heterogeneity origins. In Sec. 5, we address an important topic of first-passage time statistics that have been studied for diffusing diffusivity [38][39][40], Brownian motion subordinated by a Lévy process [41] and general random diffusivity models [39,42] but not for their combination. Based on our description, we calculate, for these mixed models, the first-passage time statistics in a bounded domain and on the real half-axis.

Integrated-diffusivity zoology

In this section we make a non-exhaustive inventory of models of non-Gaussian diffusion that can be described via the subordination approach. For the sake of clarity, we focus on one-dimensional diffusion here but most statements remain valid in any dimension (see, e.g., Sec. 5.1). In order to compute the probability density function P (x, t) of being at x at time t starting from the origin at time 0, one specifies the probability density p(x, u) of a displacement x with integrated diffusivity u and integrate it to the probability density function R(u, t) of having u in a time t such that

P (x, t) = ∞ 0 p(x, u)R(u, t)du. (1)
While the integrated diffusivity u has units of squared length, it is sometimes convenient to rescale it by diffusivity to get units of time. Whatever the convention is used, u will always appear as a dimensionless ratio u/µ u , where µ u is the scale of integrated diffusivity.

Therefore, in the present form, u can be interpreted equivalently as an integrated diffusivity, as an operational time, or as a superposition of independent processes. We will use this property in Sec. 4, to combine different models of disorder.

In the case of one-dimensional Brownian motion, the diffusion coefficient is constant so the integrated diffusivity is deterministic u = µ u t/µ t , where µ u and µ t are respectively the scales of integrated diffusivity and of physical time, from which we have

R(u, t) = δ(u -µ u t/µ t ), (2)
and the distribution of displacement reads

p(x, u) = 1 2 πµ 2 x u/µ u exp - x 2 4µ 2 x u/µ u , (3)
with the length-scale µ x . After integration in Eq. ( 1), we retrieve the Gaussian propagator

P (x, t) = 1 2 √ πDt exp - x 2 4Dt , (4)
where the diffusion coefficient D = µ 2 x /µ t is the ratio of the squared length-scale µ 2

x to the time-scale µ t . Note that the scale associated to u disappears after integration because we chose the same scale µ u in Eqs. (2,3); along the rest of the paper the final propagators will include only µ x and µ t due to the consistent choice of intermediate scales.

Time-fractional diffusion

Perhaps the best studied model of annealed disorder leading to non-Gaussian distribution is the one of time-fractional diffusion of order β. This anomalous diffusion process that displays weak ergodicity breaking can be obtained as the long-time limit of a continuous time random walk in which the waiting times are power-law distributed with an infinite mean. The probability density function of displacements is expressed in terms of the Fox H function by using a Mellin transform of the time-fractional diffusion equation [43].

The very same results can equivalently be obtained by subordination of Brownian motion with an inverse Lévy stable subordinator [34,37]. The pathwise representation to timefractional diffusion equation is

x(t) = y(S(t)), (5)
with y(t) the standard Brownian motion (or Wiener process) and S(t) the subordinator

S(t) = inf{τ > 0 : U β (τ ) > t}, (6)
where U β (τ ) is the stable subordinator with an exponent 0 < β < 1. The associated propagator P (x, t) takes the form of a subordination integral

P (x, t) = ∞ 0 p(x, τ )T (τ, t)dτ, (7)
where p(x, τ ) plays the same role as p(x, u) in Eq. (1) (i.e. the probability density of displacements) and T (τ, t) is the probability density of the integrated operational time

T (τ, t) = ∂ ∂τ P    t 0 S(t ′ )dt ′ < τ    . (8)
In the case of time-fractional diffusion, the probability density T (τ, t) of the inverse stable subordinator is given by [34] T

(τ, t) = t βτ 1+1/β µ 1/β τ µ t f β tµ 1/β τ µ t τ 1/β , (9)
where µ τ is the time-scale of internal time and f β (z) is a one-sided Lévy distribution defined by the inverse Laplace transform of fβ (s) = exp -s β . The exact form of f β (z) is known under different equivalent forms

f β (z) = 1 z H 1,0 1,1       z -β (0, β) (0, 1)       = 1 β H 1,0 1,1       z -1 (1, 1) (1/β, 1/β)       = 1 βz 2 H 1,0 1,1       z -1 (-1, 1) (-1/β, 1/β)       , (10)
where

H 1,0 1,1    x (a, A) (b, B)    is the Fox H function.
The first and second forms can be found in [START_REF] Mathai | The H-Function, Theory and applications[END_REF] and the last in [START_REF] Barkai | [END_REF] and references therein. As it simplifies Mellin transform calculations (see Sec. 4.2), the first form will be used in the rest of the paper. One gets thus

T (τ, t) = 1 βτ H 1,0 1,1    µ β t τ t β µ τ (0, β) (0, 1)    . (11)

Diffusing diffusivity

The subordination concept has also been used to describe the propagator P (x, t) of a diffusion process with time-dependent diffusivity [18,19]. In this case, P (x, t) is calculated by another subordination integral

P (x, t) = ∞ 0 p(x, u)Q(u, t)du, (12)
where Q(u, t) is interpreted as the probability density of integrated diffusivity

Q(u, τ ) = ∂ ∂u P    t 0 D(t ′ )dt ′ < u    . ( 13
)
Several models have been explored in which the mean squared displacement is linear but the diffusivity is fluctuating around its mean [14][15][16][17][18][19]. For instance, diffusivity was modeled as the distance from the origin of an n-dimensional Ornstein-Uhlenbeck process [14,[16][17][18], as a Cox-Ingersoll-Ross (CIR) process [19], as an Ornstein-Uhlenbeck process with right-sided Lévy noise [46], with stretched exponential distribution of diffusivity [20] or as the square of a Gaussian process with a non-monotonously decaying autocorrelation function [22]. The Langevin equation associated to positive time-dependent diffusivity D t can be written

dD t = a D (D t )dt + b D (D t )dW t , (14)
Y t = t 0 D t ′ dt ′ . ( 15
)
The Laplace transform of the probability density of Y t conditioned on the value

D 0 of D t at t = 0, Q(s, t|D 0 ) = E   exp   -s t 0 D t ′ dt ′   D 0   , (16)
satisfies the Feynman-Kac formula which for the considered case of time-independent coefficients a D and b D can be written via time reversal as [42] -

∂ ∂t + a D (D 0 ) ∂ ∂D 0 + 1 2 b 2 D (D 0 ) ∂ 2 ∂D 2 0 -sD 0 Q(s, t|D 0 ) = 0, (17)
subject to the condition Q(s, t = 0|D 0 ) = 1. The Laplace transform of the probability density of the integral of the CIR process Q(D, t|D 0 ) was named Υ(D, t|D 0 ) in our previous paper [39]; here we keep the notation Υ(D, t|D 0 ) for the Laplace transform of the probability density of any integrated diffusivity process that will be used in the context of first-passage time statistics in Sec. 5.1. The marginal probability density is obtained by averaging the initial diffusivity D 0 over the stationary distribution of diffusivity Π(D 0 )

Q(s, t) = ∞ 0 Π(D 0 ) Q(s, t|D 0 ) dD 0 . ( 18
)
Assuming the Gaussian form (3) of p(x, u) in Eq. ( 12) and performing the Fourier transform with respect to x, one gets

P (k, t) = ∞ -∞ e ikx P (x, t)dx = ∞ 0 e -µ 2 x k 2 u Q(u, t)du = Q(s = µ 2 x k 2 , t). ( 19
)

Superposition of independent processes

Either when the observed dynamics is sampled from a collection of independent particles with a constant but random diffusion coefficient, or when describing the particles motion in a spatially heterogeneous medium within a displacement range smaller than the correlation length of disorder, an appropriate description is obtained through a superposition of processes which in the former case could be called heterogeneous ensemble of non-Gaussian particles [25] while the latter case is termed superstatistics [18,19,26,27]. The propagator of the ensemble of particles, P (x, t), is a superposition of individual propagators p(x, t|D)

with diffusion coefficients picked from a probability density Π(D). When the spatial coordinates of the individual propagator p(x, t|D η ) scale as ∝ D η t η (i.e. anomalous diffusion) with η ∈ [0, 2] being the anomalous exponent and D η the generalized diffusion coefficient, then the propagator can be written p(x|D η t η ). The superimposed propagator,

P (x, t) = ∞ 0 p(x|D η t η )Π(D η )dD η , (20)
can be reformulated in terms of a subordination integral through the variable change

u = D η t η P (x, t) = ∞ 0 p(x, u)R(u, t)du, (21)
where the probability density of integrated diffusivity is

R(u, t) = 1 t η Π u t η . ( 22
)

Continuous-time random integrated diffusivity (CTRID)

The probability density of integrated diffusivity in the case of diffusing diffusivity presented Sec. 2.2 can be difficult to handle due to the intricate time dependence of the propagator. In order to be able to combine fluctuating diffusivity with other models of heterogeneity, we propose a simple model of integrated diffusivity that reproduces the main features of diffusing diffusivity. Instead of describing the fluctuations of diffusivity, we directly model its integral. We model the integrated diffusivity by a piecewise constant function which takes increasing random values on successive time intervals of random durations. We define the transit times to be the random time spent between two successive increments of integrated diffusivity. We consider both the transit times and the increments of integrated diffusivity to be independent, although this assumption can be relaxed (see [47]). While the piecewise constant integrated diffusivity and a continuously increasing one are different at short times, their long-time behaviors are expected to be the same. Using the renewal theory, we construct a formula for the probability density of integrated diffusivity R(u, τ ) (see Eq. ( 55)). Then, based on the long-time behavior of our version of the Montroll-Weiss formula, we establish an advection equation and a fractional advection equation that describe long-time and large-integrated-diffusivityincrements behavior of the CTRID that we solve and apply in next sections.

We also note that our CTRID model is different from the annealed-transit-time model in which the diffusivity is piecewise constant [23]. Figure 1 illustrates a random path of integrated diffusivity in diffusing diffusivity model (DD), ATTM model and our CTRID model, suggesting that all these models should be equivalent in the long-time limit. In the two latter cases, the parameters of the models can be tuned to be much closer to the DD model than presented here.

Renewal formula

We obtain the distribution R(u, τ ) by adapting the continuous-time random walk framework. We choose a stationary strictly positive distribution of integrated diffusivity jumps with density r(u) and a stationary strictly positive distribution of transit times with density φ(t). The probability density for an integrated diffusivity to be equal to u after m jumps is which in the Laplace domain becomes

r m (u) = u 0 r m-1 (u -u ′ )r(u ′ )du ′ , (23)
rm (s u ) = [r(s u )] m . (24)
In the same way, one obtains the probability density of making m integrated diffusivity increments in a time τ in the Laplace domain

φm (s τ ) = [φ(s τ )] m , (25)
whereas the probability of no integrated diffusivity change for a time τ reads in the Laplace domain as

Φ0 (s τ ) = 1 -φ(s τ ) s τ . ( 26
)
Summing over all possible numbers of jumps, one gets for the double Laplace transform of R(u, τ ):

R(s u , s τ ) = ∞ m=0 rm (s u ) φm (s τ ) Φ0 (s τ ), (27)
which simplifies to an analog of Montroll-Weiss formula

R(s u , s τ ) = 1 -φ(s τ ) s τ 1 1 -r(s u ) φ(s τ ) . ( 28
)
For instance, when r(s u ) is a Gamma distribution, one obtains the integral of a timedelayed version of the Gamma process studied in [35].

Advection equation

We first consider the case when both transit times and increments of integrated diffusivity have finite means µ τ and µ u respectively. For small s τ and s u , one has

φ(s τ ) = 1 -µ τ s τ + . . . (29) and r
(s u ) = 1 -µ u s u + . . . , (30)
so that R(s u , s τ ) ≈ µ τ µ u s u + µ τ s τ . (31)
Using the properties of the Laplace transform we obtain that in the long-transit-time and large-integrated-diffusivity-jump limit, R(u, τ ) is described by an advection equation

∂ ∂τ R(u, τ ) = - µ u µ τ ∂ ∂u R(u, τ ). ( 32
)
Given an initial distribution R 0 (u) that verifies R 0 (u) = 0 for u < 0, the solution is

R(u, τ ) = R 0 u - µ u µ τ τ . (33)
The functional form R 0 (z) of the distribution is preserved but the mode of the distribution has a constant drift over time. In this way, u is strictly increasing and it is thus an appropriate subordinator. The particular choice R 0 (z) = δ(z) makes the integrated diffusivity deterministic, implying a constant diffusion coefficient. The integrated diffusivity is characterized by both a diffusivity scale µ u representing the strength of integrated diffusivity fluctuations, and a timescale µ τ .

We illustrate the behavior of this model by considering a simple case when R 0 (u) is

an exponential distribution R 0 (z) = 1 µ u e -u/µu Θ(u), (34)
with the integrated diffusivity scale µ u , and Θ(u) is the Heaviside function. The Laplace

transform of R(u, t) with respect to u is R(s u , t) = 1 1 + µ u s u exp - µ u µ τ s u τ . (35)
Using Eq. (3) for p(x, u) we deduce the characteristic function

P (k, t) = 1 1 + µ 2 x k 2 exp - µ 2 x µ t k 2 t . ( 36
)
The first factor is the characteristic function of the Laplace distribution while the second factor is that of the Gaussian distribution. When µ 2 x k 2 ≪ 1 the distribution is well approximated by a Gaussian distribution, while when µ 2

x µt t ≪ 1, it is well approximated by a Laplace distribution. When t/µ t ≫ 1, the distribution is Gaussian with exponential tails. However, for sufficiently long time, the exponential tails cannot be observed in practice due to a limited amount of measurable data. Note that this is exactly the propagator for Brownian motion with additional Laplace distributed measurement noise.

The propagator in real space is then obtained through the convolution integral

P (x, t) = 1 2µ 2 x 2t/µ t ∞ -∞ exp - |x ′ | µ x exp - (x -x ′ ) 2 4µ 2 x t/µ t dx ′ , (37)
which yields

P (x, t) = 1 2µ x π 2 exp t µ t - x µ x erfc t µ t 1 - xµ t 2µ x t + exp 2x µ x erfc t µ t 1 + xµ t 2µ x t , (38)
where erfc(x) is the complementary error function. Figure 2 illustrates the propagator from Eq. ( 38) for CTRID with Gaussian distributed displacements. As expected, the propagator is close to a Laplace probability density function at short time and then converges to a Gaussian one a longer time. Similar calculations can be conducted for more general R 0 (u) (e.g. generalized Gamma distribution [20]).

Because the distribution of the integrated diffusivity preserves its shape at all times, the mean squared displacement, as well as all even moments of displacements, are not zero at time t = 0. For this reason, we characterize the distribution by its moments minus their initial values. In this case the mean squared displacement is a linear function of time,

X 2 (t) -X 2 (0) = 2 µ 2 x µ t t,
and the fourth moment is a quadratic function of time

X 4 (t) -X 4 (0) = 12 µ 2 x µ t t 2 + 24 µ 4 x µ t t, (40)
such that the shifted normalized excess kurtosis is

γ(t) = 1 3 X 4 (t) -X 4 (0) ( X 2 (t) -X 2 (0) ) 2 = 2 µ t t . (41)
The long-time behavior is the same as in the diffusing diffusivity model at times longer that the diffusivity autocorrelation time [19]. As this minimalistic model reproduces well the convergence to Gaussian probability density of displacements of a process with uncorrelated fluctuations of diffusivity, it can be considered as a simplified model of 'diffusing diffusivity'.

The fractional advection equation

Now we get to the case of heavy-tailed distribution of both transit times and integrated diffusivity jumps. The developments for small s τ and small s u read:

φ(s τ ) = 1 -(µ τ s τ ) δ + . . . (42) and r
(s u ) = 1 -(µ u s u ) γ + . . . , (43)
where µ τ and µ u are respective scales, while 0 < δ < 1 and 0 < γ < 1 are the scaling exponents. Eq. ( 28) yields thus

R(s u , s τ ) ≈ µ δ τ s δ-1 τ µ δ τ s δ τ + µ γ u s γ u , (44)
so that R(u, τ ) satisfies the fractional advection equation

∂ ∂τ R(u, τ ) = - µ γ u µ δ τ 0 D 1-δ τ [ 0 D γ u R(u, τ )] , (45)
where 0 D 1-δ τ and 0 D γ τ are the Riemann-Liouville operators, e.g.

0 D 1-δ τ f (τ ) = 1 Γ(δ) ∂ ∂τ τ 0 1 (τ -τ ′ ) 1-δ f (τ ′ )dτ ′ . ( 46
)
The inverse Laplace transform of Eq. ( 44) with respect to s τ yields [48] R(s u , τ

) = E δ - µ γ u µ δ τ s γ u τ δ , (47)
where E δ (z) is the Mittag-Leffler function. The inverse Laplace transform of Eq. ( 47) with respect to s u can be expressed in terms of the Fox H function [START_REF] Mathai | The H-Function, Theory and applications[END_REF] R(u, τ

) = 1 u H 1,1 2,2    u γ µ δ τ µ γ u τ δ (1, 1), (1, δ) (1, 1), (1, γ)    . (48)

Multiple subordination

In this section we discuss the possibility of describing time-dependent heterogeneities of mixed origins. The Langevin equation associated with time-dependent diffusivity and waiting times between displacements can be written as a system of three independent Langevin equations

             dx du = √ 2 dη(u) du dτ = D(τ ), dτ dt = S(t), (49)
where dη(u) is a continuous noise term (e.g. increments of a Lévy process) with a correlation function dη(u) dη(u ′ ) = δ(u-u ′ ), D(τ ) is the time-dependent diffusivity, and S(t) describes the increments of the internal time; both D(τ ) and S(t) can be stochastic processes. There are at least two methods to solve these equations. The first method is to consider the three-dimensional Fokker-Planck equation for the joint probability density P (x, u, τ, t). The marginal propagator P (x, t) that can be estimated from an experiment, is then obtained by averaging over u and τ

P (x, t) = ∞ 0 ∞ 0 P (x, u, τ, t) du dτ. (50)
We choose the second method consisting in solving three independent Fokker-Planck equations: (i) one on the probability density p(x, u) of displacements, (ii) the other on the probability density R(u, τ ) of the integrated diffusivity and (iii) the last one describing the evolution of the probability density T (τ, t) of the internal time of waiting times. The mathematical key is that a subordinated subordinator is itself a subordinator. Therefore, the probability density obtained by successively subordinating Brownian motion with different subordinators can be seen as Brownian motion subordinated by a single more sophisticated subordinator. For example, the pathwise properties of a particle with Gaussian displacements with diffusing diffusivity and waiting times between jumps can be represented as

x(t) = y[I(S(t))], (51)
with the standard Brownian motion y(t) and the subordinator I(t) interpreted as integrated diffusivity

I(t) = inf    τ > 0 : µ t µ 2 x τ 0 D s ds > t    . ( 52
)
In turn S(t) is the subordinator associated to waiting times, defined in Eq. ( 6), whose distribution is given in Eq. ( 11). Let us define x(t) = y ⋆ S 1 ⋆ S 2 ⋆ (t) as the operation associated to a double subordination

x(t) = y[S 1 (S 2 (t))].
In general, the subordination can be done an arbitrary number of times, and the characteristic function of the n-fold subordinated process

x(t) = y ⋆ S 1 ⋆ S 2 ⋆ . . . S n ⋆ (t), (53)
with y(t) distributed according to p(x, u), is

P (x, t) = ∞ 0 ∞ 0 . . . ∞ 0 p(x, u 1 )R 1 (u 1 , u 2 )R 2 (u 2 , u 3 ) . . . R n (u n , t) du 1 du 2 . . . du n (54)
with the probability density R i of the i-th subordinator.

In this section we study the propagator for processes with combined time-fluctuating diffusivity and waiting times between jumps by combining Eq. ( 7) and Eq. ( 12) into a double subordination formula

P (x, t) = ∞ 0 ∞ 0 p(x, u)R 1 (u, τ )R 2 (τ, t) du dτ. ( 55
)
This representation allows one to model a great variety of diffusion phenomena by adjusting p(x, u), R 1 (u, τ ) and R 2 (τ, t).

In Table 1, we list some examples of diffusion processes that can be described within the double subordination framework. In each case we present the characteristic function of displacement p(k, u), and the probability densities R 1 (u, τ ) and R 2 (τ, t) of two subordinators, necessary to get the characteristic function of the desired mixed non-Gaussian model. Note that the double subordination can also be used to model phenomena with two-power-law relaxation responses such as dielectric breakdown [50,51].

Superposition of processes with Continuous Time Random Integrated Diffusivity

When analyzing single particle data from an experiment on living cells, a representative statistical set is often obtained from numerous cells. If we assume that the dynamics of an individual tracer can be described by a diffusing-diffusivity model, then each tracer may have a different long-term average diffusion coefficient due to the cell-to-cell variability.

Then the propagator obtained from the ensemble of acquired trajectories is a superposition of diffusing diffusivity processes. Taking the CTRID model as a example of diffusing diffusivity process, our approach allows for the first time to describe such superposition:

P (k, t) = ∞ 0 ∞ 0 p(k, u) R 0 u - µ u τ µ τ Θ u - µ u τ µ τ 1 t Π τ t du dτ. (56)
Setting p(k, u) = exp (-k 2 µ 2 x u/µ u ), R 0 (u) to be defined in Eq. ( 34), and the distributed diffusion coefficient across the population of cells to follow a Gamma distribution

Π (z) = z ν-1 µ ν t Γ(ν)µ ν τ exp - zµ t µ τ , (57)
with dimensionless scale µ t /µ τ and shape parameter ν, the characteristic function is

P (k, t) = 1 (1 + µ 2 x k 2 ) (1 + µ 2 x k 2 t/µ t ) ν . (58)
At this stage, we notice that the initial non-Gaussian propagator converges to another non-Gaussian propagator in the long-time limit. The cross over between two non-Gaussian regimes is an important feature that only arises when at least two heterogeneity sources are present.

Model p(k, u) R 1 (u, τ ) R 2 (τ, t) Brownian motion exp -k 2 µ 2 x u/µ u δ(u -µ u τ /µ τ ) δ (τ -µ τ t/µ t ) Random diffusivity (CTRID) exp -k 2 µ 2 x u/µ u R 0 u - µ u τ µ τ Θ u - µ u τ µ τ δ (τ -µ τ t/µ t ) Superposition of Brownian motions exp -k 2 µ 2 x u/µ u 1 τ Π u τ δ (τ -µ τ t/µ t ) Diffusing diffusivity exp -k 2 µ 2 x u/µ u Q(u, τ ) δ (τ -µ τ t/µ t ) Superposition of diffusing diffusivities exp -k 2 µ 2 x u/µ u Q (u, τ ) 1 t Π τ t Space-fractional diffusion (Symmetric) exp (-|k| α µ α x u/µ u ) δ(u -µ u τ /µ τ ) δ (τ -µ τ t/µ t ) Time-fractional diffusion exp -k 2 µ 2 x u/µ u δ(u -µ u τ /µ τ ) 1 βτ H 1,0 1,1 µ β t τ t β µ τ (0, β) (0, 1) Diffusivity-fractional diffusion exp -k 2 µ 2 x u/µ u 1 u H 1,1 2,2 u γ µ δ τ µ γ u τ δ (1, 1), (1, δ) (1, 1), (1, γ) δ (τ
-µ τ t/µ t ) Space-diffusivity-time fractional diffusion exp (-|k| α µ α x u/µ u ) 1 u H 1,1 2,2 u γ µ δ τ µ γ u τ δ (1, 1), (1, δ) (1, 1), (1, γ)
1 βτ H 1,0 1,1 µ β t τ t β µ τ (0, β) (0, 1)
Table 1. Summary of models of non-Gaussian diffusion of mixed origins that can be obtained by double subordination. R 0 (u) is the distribution of integrated diffusivity obtained in Eq. ( 33), Π(z) is the distribution of diffusion coefficients distributed across a population of diffusing particles (see Eq.( 22)), Q (u, τ ) is the probability density of integrated diffusing diffusivity (see Eq.( 18)), the seventh and the last rows use for R 2 (τ, t) the probability density of the inverse stable subordinator defined in Eq (11) and the two last rows use for R 1 (u, τ ) the solution to the fractional advection equation (48).

The distribution is then the convolution of a Laplace distribution with the probability density ( 16) from [19]

F -1   1 1 + µ 2 x µt k 2 t ν   = 2 1/2-ν Γ(ν) νµ t πµ 2 x t |x| µ x νµ t t ν-1/2 K ν-1/2 |x| µ x νµ t t , (59)
where K ν (z) is the modified Bessel function of the second kind.

The propagator in real space can be computed through the convolution integral

P (x, t) = 2 1/2-ν µ 2 x Γ(ν) νµ t πt ∞ -∞ exp - |x -x ′ | µ x |x ′ | µ x νµ t t ν-1/2 × K ν-1/2 |x ′ | µ x νµ t t dx ′ , (60)
which in the case ν = 1 reduces to

P (x, t) = 1 µ x t µt -1 2 π t µ t exp - |x| µ x µ t t -exp - |x| µ x , (61)
which remains positive and properly normalized at all times. Figure 3 illustrates the propagator from Eq. ( 60), which behaves at short times as a Laplace distribution and then converges to an even less Gaussian distribution at long time (the kurtosis has increased).

Space-diffusivity-time fractional diffusion

Diffusion processes with a stable distribution of particle jumps (space fractional diffusion) and power law waiting times between jumps (time fractional diffusion) have been investigated [52][53][54][55][56]. Here our approach allows us to generalize these processes by introducing an example of non-Gaussian diffusion in which four mechanisms leading to anomalous diffusion are superimposed that we call "space-diffusivity-time fractional diffusion".

The jumps are distributed according to a symmetric stable distribution and exhibit fluctuating diffusivity for which the integral is modeled according to the fractional advection equation developed in Sec. 3.3. Finally there are power law waiting times between displacements of the particle similar to that of the time-fractional diffusion model.

The system of integro-differential equations describing the whole system is

(62a) ∂ ∂u p(x, u) = µ α x µ u ∂ α ∂|x| α p(x, u), (62b) ∂ ∂τ R(u, τ ) = - µ γ u µ δ τ 0 D 1-δ τ [ 0 D γ u R(u, τ )] (62c) 0 D β t T (τ, t) = - µ τ µ β t ∂ ∂τ T (τ, t),
where the first equation describes diffusion with symmetric stable displacements, the second one is the fractional advection equation for the integrated diffusivity, and the last one is the time-fractional advection equation of the inverse stable subordinator [34]. The characteristic function (i.e., the Fourier transform with respect to space coordinate of the solution to Eq. (62a)),

p(k, u) = exp - µ α x µ u |k| α u (63)
is first combined with the probability density of integrated diffusivity from Eq. ( 48) such that

Ṽ (k, τ ) = ∞ 0 p(k, u)R(u, τ ) du (64)
is the Laplace transform of R(u, τ ) with respect to u with parameter s u = µ α x |k| α /µ u :

Ṽ (k, τ ) = E δ - µ αγ x µ δ τ |k| αγ τ δ , (65)
that can be written in terms of the Fox H-function

Ṽ (k, τ ) = H 1,1 1,2    µ αγ x µ δ τ |k| αγ τ δ (0, 1) (0, 1), (0, δ)
   . (66)
The probability density T (τ, t) of the inverse stable subordinator is given by Eq. (11). Then the characteristic function of displacements at time t is expressed as the subordination integral

P (k, t) = 1 β ∞ 0 τ -1 H 1,1 1,2    µ αγ x µ δ τ |k| αγ τ δ (0, 1) (0, 1), (0, δ)    H 1,0 1,1    µ β t τ t β µ τ (0, β) (0, 1)    dτ, (67)
which can be interpreted as a Mellin transform ∞ 0 τ s-1 f (τ )dτ with parameter s = 0 of the product of two Fox H functions [START_REF] Mathai | The H-Function, Theory and applications[END_REF], which is

P (k, t) = 1 β H 1,2 2,3    µ αγ x µ δβ t |k| αγ t δβ (1, δ), (0, 1) (0, 1), (1, δβ), (0, δ)    . (68)
Then, the inverse Fourier transform is performed by using the properties of the Fox Hfunction which gives access to the exact expression for the probability density

P (x, t) = 1 β|x| √ π H 3,1 3,4    µ δβ t t δβ |x| 2µ x αγ
(1, 1), (0, δβ), (1, δ)

(1/2, αγ/2), (0, δ), (1, 1), (1, αγ/2)

   . (69)
Because α and γ appear as a product in the propagator, the changes of fractional integrated diffusivity plays the same role as the stable displacements of a Lévy flight. In turn, although β and δ appear as a product that scales time, they differ in the arguments of the Fox H function, suggesting new effects from the combination of both mechanisms.

The combined effect of β < 1 and δ < 1 is to slow down the dynamics while the combined effect of α < 2 and γ < 1 is to make more probable large displacements. For the first time, we are able to show that the combination of waiting times between jumps, integrated diffusivity jumps and transit-times between these, when they are all power-law distributed, leads to a scaling distribution that does not change its shape with respect to time.

First-passage times

We now look into the problem of the first-passage time statistics which is central for the description of reaction kinetics.

First-passage time in any bounded domain in R d

In the case of first-passage time in a bounded domain Ω in R d , many results are known for Brownian motion. The approach developed in this article was to design an elaborate annealed disorder description of non-Gaussian diffusion which can be shown to be equivalent to a subordinated Brownian motion. As a consequence, as explicited in [39],

the first-passage time statistics can be obtained by applying appropriate subordination integrals to the similar statistics for Brownian motion. The propagator of the process inside a bounded domain starting at a position x 0 ∈ Ω is given by

P Ω (x, t|x 0 ) = ∞ m=0 u m (x)u m (x 0 )Υ(λ m , t), (70)
where λ m is the m-th eigenvalue and u m (x) is the corresponding L 2 -normalized eigenfunction of the Laplace operator ∆, i.e. solutions of the eigenvalue problem for normal diffusion ∆u m + λ m u m = 0 in Ω subject to appropriate boundary conditions (e.g.

Dirichlet or Neumann). Here, the subscript 'Ω' emphasizes that the propagator is defined in this domain. The function Υ(λ m , t) is the Laplace transform of the probability density of integrated diffuvity in Eq. ( 54) given by the relation

Υ(λ, t) = ∞ 0 ∞ 0 . . . ∞ 0 e -λu 1 R 1 (u 1 , u 2 )R 2 (u 2 , u 3 ) . . . R n (u n , t) du 1 du 2 . . . du n , (71)
which can be written as combined subordination integral with respect to the Laplace transform of the first subordinator

Υ(λ, t) = ∞ 0 . . . ∞ 0 R1 (s u 1 = λ, u 2 )R 2 (u 2 , u 3 ) . . . R n (u n , t) du 2 . . . du n . ( 72
)
Then the probability density of the first-passage time is

ρ(t|x 0 ) = - ∞ m=0 u m (x 0 )   Ω u m (x)dx   ∂ ∂t Υ(λ m , t). (73)
The macroscopic flux is the average of this density with respect to the initial concentration of particles c(x 0 ) in the medium

J(t) = - ∞ m=0   Ω c(x 0 )u m (x 0 )dx 0     Ω u m (x)dx   ∂ ∂t Υ(λ m , t). (74)
This approach gives access to various first-passage time statistics for any combination of time-dependent heterogeneities.

For instance, in the case of Diffusivity-Time Fractional diffusion model which is the solution of Eq. ( 68) with α = 2 to ensure the process to be a subordinated Brownian motion, one gets

Υ(λ m , t) = 1 β H 1,2 2,3    µ 2γ x µ δβ t λ γ m t δβ
(1, δ), (0, 1) (0, 1), (1, δβ), (0, δ)

   . (75)
It is important to notice that the exponent γ does not affect the spectral properties of the Laplace operator.

First-passage times on a positive half-axis

The above analysis is not applicable to first-passage times in unbounded domains for which the Laplace operator has a continuous spectrum. In some cases, one can employ symmetries to get explicit results (see [39] for details). Here, we illustrate this possibility for diffusion on a half-line. From the characteristic function of displacements, the survival

The asymptotic behavior of the Fox H function has been investigated in [START_REF] Mathai | The H-Function, Theory and applications[END_REF]. The probability density of the first-passage times is then obtained as ρ(t|x 0 ) = -∂S(t|x 0 )/∂t.

Conclusion

In this article, we investigated the probability density P (x, t) of finding a particle, started at time 0 from the origin, at point x at time t subject to non-Gaussian diffusion of mixed origins. We showed that various physical mechanisms that can be described by subordinated Brownian motion, can be combined and investigated as a whole. The necessary step is to calculate the probability density of the integrated diffusivity (or equivalently internal time). Therefore we proposed a model of continuoustime random integrated diffusivity for which the asymptotic behavior in two regimes described by an ordinary or fractional advection equation, is explicitly derived. From the combined subordination technique and the CTRID model, we derived the propagator P (x, t) of the superposition of Gaussian processes with integrated diffusivity modeled as a CTRID as well as the space-diffusivity-time fractional model. We showed that combined heterogeneity delivers rich phenomenology; for example, the probability density of displacements can be non-Gaussian at short-time and converge to a second non-Gaussian distribution at long-times. Finally we showed how the first-passage time statistics in any bounded domain in R d (as well as in the half-line) can be obtained for the mixed models if the results are known for Brownian motion. The present work aims at shedding light on diffusion in disordered media in which several disorder mechanisms are simultaneously present.

  with some drift a D and volatility b D . The positivity of D t can be ensured either by choosing the specific form of the drift and volatility functions in the Langevin equation, or by introducing a reflecting boundary at 0 to the associated Fokker-Planck equation. To describe the fluctuating diffusivity by subordination, one needs to obtain the probability density of the integrated diffusivity

Figure 1 .

 1 Figure 1. Illustration of the continuous-time random integrated diffusivity (red dashed line) where the integrated diffusivity takes piecewise constant values u 1 , u 2 , . . . , u n during random times τ 1 , τ 2 , . . . , τ n . The black dash-dotted line shows the annealed transit time model where the diffusion coefficient D 1 , D 2 , . . . , D n remains constant during random times τ 1 , τ 2 , . . . , τ n . Finally the blue line presents the continuously increasing integrated diffusing diffusivity as a function of time (arbitrary units).

Figure 2 .

 2 Figure 2. Propagator from Eq. (38) for a CTRID model with Gaussian displacements with parameters µ x = 1 and µ t = 1 at time t = 0.01 (red line), t = 0.5 (blue dashed line) and t = 10 (black dashed-dotted line).

Figure 3 .

 3 Figure 3. The propagator of a superposition of processes with Gaussian displacements and integrated diffusivity modeled by CTRID obtained from Eq. (60) with parameters µ x = 3, µ t = 1 and ν = 0.3 at time t = 0.1 (red line), t = 3 (blue dashed line) and t = 10 (black dashed-dotted line).

  probability of a particle in [0, ∞] with absorbing boundary located at x = 0 is 0 ) P (k, t)dk.(76)When P (k, t) is obtained as a process with successive subordinations with stable symmetric displacement (see Eq. (63)), the survival probability is u, τ )T (τ, t)du dτ. (78)When R(u, τ ) is the solution (48) of the fractional advection equation obtained in Sec.3.3 and T (τ, t) is the inverse stable subordinator from Eq. (11), the calculation of the survival probability consists in two consecutive Mellin transforms of the product of two 1), (0, βδ), (1, δ) (0, δ), (1, 1), (1/2, αγ/2), (0, αγ/2)


Liubov Tupikina
email: liubov.tupikina@cri-paris.org

Denis S Grebenkov
email: denis.grebenkov@polytechnique.edu

Structural and temporal heterogeneities on networks

Keywords: Continuous time random walk, stochastic processes on graphs, diffusion, complex networks, heterogeneities, first passage time, network measures

A heterogeneous continuous time random walk is an analytical formalism for studying and modeling diffusion processes in heterogeneous structures on microscopic and macroscopic scales. In this paper we study both analytically and numerically the effects of spatio-temporal heterogeneities onto the diffusive dynamics on different types of networks. We investigate how the distribution of the first passage time is affected by the global topological network properties and heterogeneities in the distributions of the travel times. In particular, we analyze transport properties of random networks and define network measures based on the first-passage characteristics. The heterogeneous continuous time random walk framework has potential applications in biology, social and urban science, search of optimal transport properties, analysis of the effects of heterogeneities or bursts in transportation networks.

Introduction

Dynamical properties of random walks on networks are related to many vital questions, such as optimality and efficiency of road systems, internet search strategy, functioning of metabolic networks [START_REF] Montroll | [END_REF][START_REF] Barthelemy | Spatial networks[END_REF][START_REF] Cohen | [END_REF]. Dynamics on various types of random networks (graphs) [START_REF] Sood | First-passage properties of the Erdos-Renyi random graph[END_REF] have been extensively studied within the last decade and applied to describe epidemics and social processes [START_REF] Brockmann | The scaling laws of human travel[END_REF]. Many of large networks encountered in our everyday life have the so-called small-world property meaning that all nodes are closely connected [START_REF] Agliari | Random walks on deterministic scale-free networks: Exact results[END_REF][START_REF] Julaiti | Eigenvalues of normalized Laplacian matrices of fractal trees and dendrimers: analytical results and applications[END_REF][START_REF] Albert | Statistical Mechanics of Complex Networks[END_REF]. Moreover, it has been found that topological, geometric, and hydrological characteristics of a network system are directly linked to the characteristics of random walks on this network [START_REF] Redner | A Guide to First-Passage Processes[END_REF][START_REF] Havlin | Diffusion in disordered media[END_REF], such as the first passage time (FPT), i.e. the time of the first arrival of the random walk to a target node [START_REF] Bollt | What is special about diffusion on scale-free nets?[END_REF]. The first passage time characteristics are related to trapping problems, which play an important role in the control theory. Numerous problems on dynamical processes on networks [START_REF] Hwang | First Passage Time for Random Walks in Heterogeneous Networks[END_REF] include finding an appropriate placing of the trapping site in a network in order to obtain the needed trapping efficiency, search strategies etc. [START_REF] Lambiotte | Diffusion on networked systems is a question of time or structure[END_REF]. For homogeneous networks these problems have been studied extensively [START_REF] Maier | Cover time for random walks on arbitrary complex networks[END_REF][START_REF] Hudghes | Random walks in random environments[END_REF], for example, using the underlying backward equation approach [START_REF] Bollt | What is special about diffusion on scale-free nets?[END_REF]. However, the analytical approaches, describing random walks on feature-rich (when nodes or links are endowed with some attributes [START_REF] Monti | Estimating latent feature-feature interactions in large feature-rich graphs[END_REF]) heterogeneous networks, are lacking [START_REF] Bollt | What is special about diffusion on scale-free nets?[END_REF]. It has been demonstrated that in continuous domains the mean FPT can differ from the most probable FPT [START_REF] Godec | First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time[END_REF][START_REF] Grebenkov | Towards a full quantitative description of singlemolecule reaction kinetics in biological cells[END_REF], in particular, moments of FPT do not fully reflect the transport efficiency. Moreover, since real-world networks are highly heterogeneous, reliable methods for analyzing transport properties on them are required [START_REF] Valdano | Epidemic Threshold in Continuous-Time Evolving Networks[END_REF]. To our knowledge, the complete description of the first passage quantities on various networks with structural or temporal heterogeneities has not been achieved. There is a particular reason for studying FPT properties on random networks. Some random networks have shown clear advantages over other network topologies (such as complete graphs) in terms of transport optimality and scaling efficiency, quantities which characterize the growth of the FPT density with the increase of network size and network diameter [START_REF] Bollt | What is special about diffusion on scale-free nets?[END_REF]. Recently dynamics on random networks has been studied using the spectral and mean-field theories [START_REF] Grabow | Small-world network spectra in mean-field theory[END_REF][START_REF] Mieghem | Graph Spectra for Complex Networks[END_REF], which are not always applicable to heterogeneous networks.

In this article we consider Heterogeneous Continuous Time Random Walk (HC-TRW) model, in which heterogeneity is understood from both structural and temporal perspectives [START_REF] Grebenkov | Heterogeneous continuous-time random walks[END_REF]. We investigate HCTRW on various graphs to see the interplay between the dynamical random walk characteristics and structural properties of regular and random networks. The structure of the paper is the following. After presenting the model in subsection 1.1, we investigate in Section 2 the first passage time of HCTRW on regular and irregular networks with structural and temporal heterogeneities (or perturbations). In Sections 3 and 4 we discuss the results and conclude with the main findings.

Model

In HCTRW model [START_REF] Grebenkov | Heterogeneous continuous-time random walks[END_REF] a random walk moves on a graph G in continuous time, jumping from one node to another with a transition (stochastic) matrix Q whose element Q xx is the probability of jumping from the node x to x via link e xx and the travel time needed to move from x to x is a random variable drawn from the probability density ψ xx (t). The propagator of HCTRW in the Laplace domain was obtained as

Px 0 x (s) = 1 -x Qxx (s) s [(I -Q(s)) -1 ] x 0 x , (1)
where Q(s) is the Laplace transform of Q(t), the generalized transition matrix with elements Q xx (t), which are equal Q xx ψ xx (t) [START_REF] Grebenkov | Heterogeneous continuous-time random walks[END_REF]. From this formula other important quantities of the process can be derived. In particular, the probability density ρ x 0 (t) of the first passage time to a single absorbing node x a can be obtained using the renewal approach

ρx 0 (s) = Px 0 x a (s) Px a x a (s) . (2)
When there are many absorbing nodes x a , Eq. ( 2) is not applicable but there are other ways to compute FPT density.

For each considered network (see below) we construct the transition matrix Q and the generalized matrix Q(t), given the set of the travel time distributions ψ xx (t) for each link e xx . Then we perform the numerical computation of the inverse Laplace transform of ρx 0 (s) by Talbot algorithm to get the FPT density ρ x 0 (t) in time domain [START_REF] Talbot | The accurate numerical inversion of Laplace transforms[END_REF].

Results

We investigate the behavior of the FPT probability density on several networks. In each of these networks, we will select an absorbing node x a and calculate the probability density ρ x 0 (t) of the FPT to this node from a prescribed starting point x 0 . We will focus on the effects of heterogeneity by introducing a trap node x h , at which the travel time distribution is different from the remaining nodes. Although the HCTRW model allows for arbitrary travel time distributions, we restrict our analysis to the most common exponential distribution, for which ψxx (s) = 1/(1 + sτ ), where τ is the mean travel time, ψxx (s) Laplace transform of ψ xx (t). When a random walk arrives at the trap node, it is kept there for much longer times than τ . This mechanism can be realized by using the exponential distribution with a much larger mean travel time. However, to highlight to effect of the trap node (heterogeneity) x h , we will consider the distribution with infinite mean travel time, typically of the form ψx h x (s) = ψx h (s) = 1/(1 + s α τ α ), with an exponent α between 0 and 1. We will investigate how the presence of such a trap node and of eventual links avoiding this trap affect the FPT density. We start with the HCTRW on regular structures, such as one-dimensional lattice chain and finite fractals, then we study FPT properties on random networks.

HCTRW on regular graphs

We use the following setup of the HCTRW model on regular loopless graphs. In a given graph we also induce an avoiding link e x h +1x h -1 , which connects the nodes x h + 1 and x h -1 around the trap node x h . Then we compute FPT densities for the cases with and without the avoiding link.

First, we consider the HCTRW on a one-dimensional chain with one absorbing node, one reflecting node and two trap nodes either to the left or to the right from x 0 . The corresponding FPT densities (see Fig. 1) show the same long-time behavior t -1-α , discussed in [START_REF] Grebenkov | Heterogeneous continuous-time random walks[END_REF]. This result can be understood from the fact that a random walk on a one-dimensional chain with or without avoiding link will eventually get into a trap node even though an avoiding link potentially prevents a random walk from getting into a trap. As an example of regular loopless structures, we consider a generalized Vicsek fractal, G V gf , which resembles the dendrimers construction [START_REF] Blumen | Generalized Vicsek Fractals: Regular Hyperbranched Polymers[END_REF][START_REF] Liu | Laplacian spectra of a class of small-world networks and their applications[END_REF]. The fractal G V gf is constructed iteratively in a deterministic way by going from generation g to generation g + 1 with coordination number f . A Vicsek fractal graph for f = 3 and g = 3 is shown in Fig. 2 (top).

We use the HCTRW setup, where all travel time distributions are ψ(s) = 1/(1 + sτ ), except x h , at which ψx h (s) = 1/(1 + s α τ α ) with α = 0.5. We compute the FPT density on a Vicsek fractal for two cases: when the trap is placed on the shortest path between x 0 and x a and outside this path. As we observe from Fig. 2 (bottom) the short-time regime of the FPT density with avoiding link differs from that without avoiding link, whereas these cases give the same powerlaw scaling t -2d s -1-α in the long-time regime, where d s is the spectral dimension of the Viscek fractal. At the same time, it is known that the spectral properties of graph Laplacians are reflected in dynamical properties of random walks on them. In the case of the Vicsek fractal, the graph Laplacian obeys the simple scaling, determined by the spectral dimension of the Vicsek fractal: d s = ln(f +1) ln(3f +3) ≈ 0.557.

HCTRW on random networks

Here we consider two classes of random networks: Scale-Free (SF) and Watts-Strogatz (WS) model. We construct SF networks using the preferential attachment Barabasi-Albert model G(N, m, m 0 ), where N is the number of nodes, m is the number of initially placed nodes in a network and m 0 is the number of nodes, a newly added node is connected to [START_REF] Albert | Statistical Mechanics of Complex Networks[END_REF]. The preferential attachment mechanism drives the network degree distribution to obey a power law decay with the exponent γ ≈ 3 [START_REF] Der Hofstadt | Random Graphs and Complex Networks[END_REF]. The small-world Watts-Strogatz network model is denoted by

G(k, β),
where k is the average network degree and β is the rewiring parameter [START_REF] Watts | Collective dynamics of 'small-world[END_REF]. WS model is constructed from a regular ring lattice, a graph where each nodes is connected to k neighbors, k/2 on each side, where each link is then rewired with the probability β.

We numerically compute the FPT density on SF networks. In this section we demonstrate the FPT densities on one random realization of SF network, while we (Bottom) FPT densities on this Vicsek fractal with ψ(s) = 1/(sτ + 1), τ = 1 in all nodes except x h , at which ψ(s) = 1/(s α τ α + 1), α = 0.5. The trap x h is placed either on the shortest path between x 0 and xa (x h = 14), or outside (x h = 54). We also consider the case with an additional link avoiding the trap x h = 14, which corresponds to the transition matrix with a local perturbation.

also checked that results are qualitatively the same when considering an ensemble of random networks. First, we investigate the influence of the individually placed traps. For this we randomly fix x 0 , x a and place x h into two different communities of the network, i.e. groups of nodes, which are interlinked with each other more than with other nodes [START_REF] Barthelemy | Spatial networks[END_REF]: (a) a trap x h is placed on a path (maybe not the shortest) between x 0 and x a , denoted as x 0 < x h < x a , hence a random walk cannot avoid getting into the trap x h ; (b) x h is placed in another community in respect to x 0 node, such that a random walk is able to reach x a without passing through x h . Figure 3 illustrates these cases. The FPT densities for these two cases are compared with the FPT density for the homogeneous HCTRW, when all travel time distributions are exponential with the same parameter τ : ψ(s) = 1/(1 + sτ ). Note that the case (b) is the homogeneous case in respect to a node x a , since the trap node x h is not reached by a random walk. At the same time, placing a trap node in the same community with x 0 and x a changes the long-time behavior of the FPT density in comparison to the homogeneous case (Fig. 3 (bottom)). In contrast, for the case (b) the long-term behavior is the same as for the homogeneous case. We note here that, communities structure in larger networks can be less obvious than in the scheme in Fig. 3, and the HCTRW framework with placing trap nodes in different parts of a network can be used for analysis of communities and trapping efficiency.

Short-time regimes of FPT density on scale-free networks

From our numerical results for loopless network (subsection 2.1) we observe that the FPT densities can differ from each other in the short-time regime. We numerically calculate the FPT density for the HCTRW separately in two cases of SF networks: for SF with locally tree-like structure (m = 1 parameter of SF model) and with loops (m > 1). First we test the influence of the global topological properties on the FPT density and plot it for SF networks with fixed m 0 and various m parameters and homogeneous exponential travel times. Figure 4 shows that the case for m = 1 affects mostly the short-time and intermediate-time regimes, while cases with m > 1 slightly differ in the short-time regimes.

Then we fix all travel times to be independent identical exponentially distributed ψ(s) = 1/(sτ + 1), τ = 1, except in x h , at which ψx h (s) = 1/(s α τ α + 1) ν , α ∈ (0, 1), ν = 2, where an auxiliary parameter ν allows one to scatter more ρ x 0 (t) for different x 0 . The FPT densities are computed for different starting points x 0 and fixed x a on a SF network (Figs. [START_REF] Blumen | Generalized Vicsek Fractals: Regular Hyperbranched Polymers[END_REF][START_REF] Barrat | On the properties of small-world network models[END_REF]. Knowing that there are different shortest path distances between x 0 and x a , we estimate numerically the lengths of all possible shortest lengths |x 0 , x a | for fixed x a using using Dijkstra's algorithm. We further explore the relation between the short-time regime of the FPT density and the intrinsic network metrics, defined by the shortest paths between nodes. For SF network with loops there exist several shortest paths between randomly taken nodes. For example, in a SF network with m = 5, m 0 = 6, N = 100 three different lengths of the shortest path exist. Correspondingly, we clearly observe three different early time behaviors for the FPT density, while for large time the FPT densities behave similarly for different x 0 . This illustrates the fact that the longtime behavior of the FPT density does not capture some important information about the system.

Next we compare the FPT densities of SF tree and non-tree networks for different starting positions x 0 (Fig. 5 and Fig. 6). We observe that changing x 0 in a scale-free tree (Fig. 5) generally affects the whole FPT density, while in a scale-free network with m = 5 only the short-time regime is affected. This illustrates the fragility of the transport properties on tree-like structures. One of simple explanations is that for the same number of nodes N a tree network has larger variety of distances between nodes than in the case of a network with loops. Both, scale-free networks with and without loops obey the small-world effect. Further we study the effects of the small-world property on the FPT density.

FPT density on Watts-Strogatz network

Another example of a network with small-world property is the WS network model [START_REF] Watts | Collective dynamics of 'small-world[END_REF]. The average shortest path length in the WS model gives an estimate of the small-world property. For β = 0 (k-circular graph) the average shortest path length is N/4k [START_REF] Barrat | On the properties of small-world network models[END_REF][START_REF] Newman | Mean-field solution of the small-world network model[END_REF], while for β = 0 the average shortest path length is smaller than N/4k and for b = 0.2, k = 8 gives M ≈ 4, Fig. 7 (top). The FPT density for a WS network with various initial positions x 0 and fixed x a are shown in Fig. 7 (bottom). Here we consider the WS network with the rewiring probability β = 0.2, the average degree k = 8 and N = 100 nodes. The first observation about the short-time regime is that we clearly see m distinct groups, where m is the number of different shortest paths lengths in a network. Another observation is absence of a plateau in the intermediate regime, Fig. 7. For SF networks with m = 5 the intermediate regime has a characteristic plateau for t ∈ [10 1 , 10 3 ], Fig. 6, while for WS networks the plateau is not found. In section 3 we discuss in more details how one can compare FPT densities for different networks.

Analytical insights onto the FPT density

Here we discuss the analytical calculations for the short-time and long-time regimes of the FPT density for a particular case ofs HCTRW model. Let us consider a path of a random walk from x 0 to an absorbing node x a on a network with transition matrix Q. If all travel time distributions are exponentials with the mean τ , then the probability of making the path of length n is:

P n (t) = t n-1 τ n Γ (n) e -t/τ . (3)
In the short-time regime the function t n-1 dominates, while in the long-time regime the influence of the exponential e -t/τ is much more profound. The total probability of arriving at x a from x 0 at time t is:

P x 0 x a (t) = ∞ n=1 (Q n ) x 0 x abs P n (t) = e -t/τ τ ∞ n=1 (Q n ) x 0 x a (t/τ ) n-1 (n -1)! . ( 4
)
The calculations of the propagator for a tree-like network can be simplified using the property that the n th matrix power in the infinite series becomes zero matrix for n ≥ d, where d is the diameter of a network [START_REF] Petit | Random walk on temporal networks with lasting edges[END_REF]. Variations of the shortest path |x 0 -x a | between the starting point x 0 with the fixed absorbing node x a affect (Q n ) x 0 x a . The long-time behavior in finite graphs, regardless of a network topology, has generic exponential tail [START_REF] Bollt | What is special about diffusion on scale-free nets?[END_REF]. However, this result is valid for the long-term behavior on structures without traps. Let us consider a HCTRW model with induced structural or temporal perturbations, in which all nodes have a finite-moment distribution ψ(t) except one node x h with a heavy-tail distribution ψ x h (t). Then the spectral properties of the generalized transition matrix Q(t) are affected by structural (Q → Q pert ) and distributional (ψ(t) for x h → ψ x h (t)) perturbations [START_REF] Grebenkov | Heterogeneous continuous-time random walks[END_REF]. In the case of distributional perturbations the graph is kept fixed with the same transition matrix Q but the travel time distributions ψ xx (t) are edge dependent and the generalized transition matrix Q(t) changes. In the case of structural perturbations, adding or removing links affects the matrix Q itself, as the result, the stationary distribution changes. Both structural and distributional perturbations can be analyzed from the perspectives of temporal networks. In particular, one x 0 <x h <x a

x 0 <x a <x h Fig. 3 FPT densities for HCTRW on SF network G(N, m, m 0 ) for N = 100, m = 5, m 0 = 2 with two different positions of traps: x h is placed on a path x 0 < x h < xa (circles) and in another community of a network in respect to the starting node x 0 , x 0 < xa < x h (stars).

Travel time distributions are ψ(s) = 1/(sτ + 1), τ = 1, in a trap ψx h (s) = 1/(s α τ α + 1), α = 0.1. x 0 and xa are fixed.

can consider the behavior of a random walk in the HCTRW model as a simple random walk moving on continuously changing temporal network. In this temporal networks links are available in a certain time intervals, distributed according to some the probability distribution e xx (t) for an edge e xx . Here we refer to the frameworks presented in [START_REF] Ahmad | Tie-decay temporal networks in continuous time and eigenvector-based centralities[END_REF][START_REF] Valdano | Epidemic Threshold in Continuous-Time Evolving Networks[END_REF][START_REF] Lambiotte | Burstiness and spreading on temporal networks[END_REF], where continuously evolving networks are studied. Fig. 4 FPT density for SF networks with parameters m 0 = 11 and m = 1 (tree SF network), m ∈ [2, 9] (non-tree SF network) with N = 100 nodes. Travel time distribution for all the nodes is ψ(s) = 1/(sτ + 1), τ = 1. x 0 and xa are randomly chosen in each SF network.

Applications to real-world networks

Let us give an illustrative example of the HCTRW applications. If a train station with a small number of connections was shut down due to some technical problems, it would not significantly increase the number of additional trips that passengers would have to take to arrive at their final destinations. However, if a random perturbation affects a station with many connections or a connecting point between several clusters of stations, the average path length can increase dramatically. In real transportation systems the road occupancy is varying over time, as such, the link occupancy between two places x, x can be modeled by a random variable distributed with ψ xx (t). Then the perturbation in the transportation liquidity (bursts) at some node can be modeled in the HCTRW as the trap x h with heavy-tail edge travel times. Links avoiding a trap correspond to additional transport connections, e.g. with buses etc. We consider a spreading entity (infection, rumour, flux of passengers, etc.), which can be modelled by a random walk behavior, although in the case of epidemics spreading the spreading quantity is not necessarily conserved. We choose to analyze the London metro [START_REF]Wolfram mathematica database[END_REF], which shares some common features with other metro systems in the world, such as the average node degree of around 2.5 inside the core of metro network [START_REF] Roth | Evolution of subway networks[END_REF]. The London metro on Fig. 8 (top) has N = 299 stations (nowadays it has around 353 stations and nearly 400 edges). The core of the London metro network, neglecting the radial stations outside the circle, exhibits the small-world property. We introduce a trap (perturbation) to London metro to see how this affects the transportation properties on the whole network. In Fig. 8 we plot the FPT densities starting from different stations of the London metro and finishing at Piccadilly circus metro station (popular destination station). Comparing Fig. 8 (bottom) with the FPT density for SF and WS networks (Figs. 5, 6, 7) helps to get some insights. The whole metro network of London is separated into several groups of nodes according to the destination, in each group the short-time regimes of ρ x 0 (t) are similar. The long time regime is independent from the starting points (Fig. 8).

Discussions

The properties of random walks on various networks have been studied extensively [START_REF] Klafter | First steps in random walks[END_REF][START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF][START_REF] Steinbock | Distribution of shortest path lengths in a class of node duplication network models[END_REF]. Quite often the asymptotic analysis and the averaged characteristics of diffusive transport and random walk dynamics, such as mean return time and mean first-passage time, are calculated [START_REF] Bollt | What is special about diffusion on scale-free nets?[END_REF][START_REF] Maier | Cover time for random walks on arbitrary complex networks[END_REF]. Here we focused on studying the probability distribution of the first-passage time for random walk on complex networks with heterogeneities. We identified three regimes of the FPT density (Fig. 9): short-time, intermediatetime and long-time regimes. This scheme also demonstrates the observation about how the FPT density is affected by different types of perturbations. In particular, local structural perturbations (such as an inclusion of links avoiding x h ) mainly affect the short-time regime of the FPT density. The global structural perturba-tions (local perturbations of the network structure are accumulated and affect the network topology globally) generally have an effect on the intermediate regime. Finally, the distributional perturbations (when the travel times ψ xx (t) are linkdependent) are mainly affecting the long-time regimes of the FPT density. One can also consider the interplay of different types of perturbations but this goes beyond the scope of this paper. We start with the least studied short-time regime [START_REF] Bollt | What is special about diffusion on scale-free nets?[END_REF].

The short-time regime is strongly affected by the geometric network properties, such as the distance between x 0 and x a , which can be calculated, as the minimal power of the adjacency matrix of a network with nonzero element A x 0 x a : A L-1

x 0 x a = 0, A L x 0 x a = 0. As we saw from the numerical computations of ρ x 0 (t) (subsection 2.1), a local perturbation induced by an avoiding link strongly affects the shorttime regime (Fig. 2 (bottom)). In other words, local changes of a network structure affect the distribution of the shortest paths, which in turn affects the left tail of the FPT density (Figs. 5,6). In the same vain, the FPT properties are strongly affected by the metrics of a graph and not just by the node degrees, in the contrary to the common belief that general transport properties of networks depend exclusively on node degrees. At the same time, the non-linear network measure distribution is related to dynamical properties in a rather complex way and should be analyzed separately.

The intermediate-time regime of the FPT density mainly depends on the global topological features of a network (i.e. small-world property, loopless structure). For some networks the intermediate regime is less pronounced than for others (compare Figs. 5 and6). For instance, for WS model (Fig. 7), the intermediate regime between the most probable first-passage time, t mp , and beginning of a decrease of the FPT density is not so well pronounced. For the networks without traps changing the shortest path between x 0 and x a shifts the maximum of the FPT. This resembles the properties of diffusion in continuous domains [START_REF] Godec | First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time[END_REF], where t mp is proportional to a distance between the starting point and the target. Another important characteristics of the intermediate regime is the presence of a plateau, like the intermediate regime at t ∈ [10 1 , 10 3 ] for SF networks with m = 5, m 0 = 6 (Fig. 6), which is a general feature present also for FPT in continuous domains [START_REF] Grebenkov | Towards a full quantitative description of singlemolecule reaction kinetics in biological cells[END_REF].

Finally, the long-time regime is largely influenced by temporal heterogeneities, i.e. the trap nodes x h with a heavy tail distribution of travel time ψ x h x . In particular, it was identified that for simple random walks MFPT (mean first-passage time) depicts some long-time properties of the dynamics. However, MFPT neglects some important information about the process, significantly overestimating the scales of the most-probable first-passage time [START_REF] Godec | First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time[END_REF][START_REF] Grebenkov | Towards a full quantitative description of singlemolecule reaction kinetics in biological cells[END_REF]. At the same time it is known that for dynamics of a particle diffusing in a continuous domain the MFPT is of order L 2 /D, where L is the domain size, D is the diffusion coefficient [START_REF] Singer | Narrow Escape, Part II: The Circular Disk[END_REF]. In the continuous case the most probable FPT is shown to strongly depend on the starting position and is almost independent of the target properties [START_REF] Grebenkov | Towards a full quantitative description of singlemolecule reaction kinetics in biological cells[END_REF]. How the dynamical properties are changing when distributional perturbation at x h is introduced? The influence of x h on the long-time regime is illustrated on Fig. 3, where placing x h in different communities of a network changes the right tail of the FPT density. Note that the trap does not affect the geometric length of the path between x 0 and x a , but may affect the shortest-in-time path between nodes. Hence even if there are several paths between nodes there is a non-zero probability for a random walk to reach x h at long times, which then may affect the right tail of the distribution (Fig. 1). To complete our discussions on this point, we recall that the properties in the long-time regime of random walks on networks without heterogeneities have been proven to depend on the dimension of a random walk [START_REF] O'shaugnessy | Analytical Solutions for Diffusion on Fractal Objects[END_REF]. Moreover, in the case of the compact random walk exploration, e.g. HCTRW on Viscek fractals with spectral dimension d s < 2 (Fig. 2), the asymptotic prop- erties of the FPT density do not depend on a degree of a target node and have the power law decay t -2d s -1-α .

Conclusions

In this article we considered the HCTRW framework for studying the featurerich networks, in particular, we studied effects of heterogeneity on random walk dynamics on several types of networks. We explored FPT densities on regular and irregular structures, such as regular fractals, random Scale-Free and Watts-Strogatz network models. The comparison of regular versus irregular structures allowed us to refine some properties of the FPT densities on complex networks. In particular we analyzed effects of structural and distributional heterogeneities using the first-passage time as one of the key indicators of how fast information diffuses in a given system [START_REF] Maier | Cover time for random walks on arbitrary complex networks[END_REF][START_REF] Redner | A Guide to First-Passage Processes[END_REF]. Heterogeneities are encoded in the generalized transition matrix Q(t) [START_REF] Grebenkov | Heterogeneous continuous-time random walks[END_REF], which affects its spectral properties and as the result the dynamical properties of processes on such networks. The influence of the structural and distributional perturbations of network structure on FPT density is summarized in Fig. 9. The local and global structural perturbations mainly influence the short-time and intermediate-time regimes. The distributional perturbations generally affect the long-time regime. Notably, not only the topology of the network but also temporal discrepancy encoded in the distributional heterogeneities ψ xx (t) can significantly alter the behavior of random walk on a network.

In particular, we found that the short-time regime depends on the distribution of distances between starting and absorbing nodes in a network. The existence of loops generally influences the dynamical properties of random walks on these graphs [START_REF] Dorogovtsev | Evolution of networks[END_REF]. For the HCTRW model the changes of global network structure influence the intermediate-time regime of the FPT density, e.g. SF network with or without cycles has different intermediate regime duration (section 2). In general, analysis of the short-time and intermediate-time regimes allow us to characterize transport properties of networks even when the first moment of the distribution is infinite [START_REF] Hernandez-Garcia | First-passage-time statistics in disordered media[END_REF].

To summarize, the FPT density provides us with alternative characteristics of network topology in addition to known static measures, such as degree or closeness network measures [START_REF] Barthelemy | Spatial networks[END_REF]. The average number of distinct sites, visited by HCTRW and other observables derived from the FPT densities can be a quantity with the practical relevance for charaterization of spreading processes.

The presented analysis of random walk dynamics can be further explored on real-world networks, as it has been done with the CTRW approach for analysis of human travel laws in the continuous domains [START_REF] Brockmann | The scaling laws of human travel[END_REF]. We applied the HCTRW model on a network graph of the London metro. The next step would be to compare the results from HCTRW simulations, Fig. 8, with the original data of travelling passengers on real-world networks [41]. The transportation network is the network type, to which the HCTRW framework can be further applied, taking into account their nontrivial stochastic nature. Moreover, the HCTRW framework can be used

FIGURE LEGEND

Figures:

Fig. 1 FPT density on a chain graph with N = 100 nodes, ψ(s) = 1/(sτ +1), τ = 1 in all nodes except the trap nodes either at x h = 25, or at x h = 75, at which ψx h (s) = 1/(s α τ α + 1), α = 0.5, x 0 = 50 and x a = 1. The additional link is placed around the trap x h = 25. Fig. 2 (Top) Vicsek fractal G V 33 with N = 64 nodes. The starting node x 0 is in the center (red circle), the absorbing node x a is one of the dead-ends (black circle), traps x h (empty circles). (Bottom) FPT densities on this Vicsek fractal with ψ(s) = 1/(sτ + 1), τ = 1 in all nodes except x h , at which ψ(s) = 1/(s α τ α + 1), α = 0.5. The trap x h is placed either on the shortest path between x 0 and x a (x h = 14), or outside (x h = 54). We also consider the case with an additional link avoiding the trap x h = 14, which corresponds to the transition matrix with a local perturbation. Fig. 3 FPT densities for HCTRW on SF network G(N, m, m 0 ) for N = 100, m = 5, m 0 = 2 with two different positions of traps: x h is placed on a path x 0 < x h < x a (circles) and in another community of a network in respect to the starting node x 0 , x 0 < x a < x h (stars). Travel time distributions are ψ(s) = 1/(sτ + 1), τ = 1, in a trap ψx h (s) = 1/(s α τ α + 1), α = 0.1. x 0 and x a are fixed. Fig. 4 FPT density for SF networks with parameters m 0 = 11 and m = 1 (tree SF network), m ∈ [2, 9] (non-tree SF network) with N = 100 nodes. Travel time distribution for all the nodes is ψ(s) = 1/(sτ + 1), τ = 1. x 0 and x a are randomly chosen in each SF network.

Fig. 1 FPT

 1 Fig.1FPT density on a chain graph with N = 100 nodes, ψ(s) = 1/(sτ +1), τ = 1 in all nodes except the trap nodes either at x h = 25, or at x h = 75, at which ψx h (s) = 1/(s α τ α + 1), α = 0.5, x 0 = 50 and xa = 1. The additional link is placed around the trap x h = 25.

Fig. 2 (

 2 Fig.2(Top) Vicsek fractal G V 33 with N = 64 nodes. The starting node x 0 is in the center (red circle), the absorbing node xa is one of the dead-ends (black circle), traps x h (empty circles). (Bottom) FPT densities on this Vicsek fractal with ψ(s) = 1/(sτ + 1), τ = 1 in all nodes except x h , at which ψ(s) = 1/(s α τ α + 1), α = 0.5. The trap x h is placed either on the shortest path between x 0 and xa (x h = 14), or outside (x h = 54). We also consider the case with an additional link avoiding the trap x h = 14, which corresponds to the transition matrix with a local perturbation.

Fig. 5 (

 5 Fig. 5 (Top) Example of SF network with a fixed xa marked in black circle. (Bottom) FPT densities on SF network G(N, m, m 0 ) N = 100, m = 1, m 0 = 6 for starting points x 0 and fixed xa such that: |x 0 -xa| = 1 (black), |x 0 -xa| = 2 (blue), |x 0 -xa| = 3 (red) (FPT densities with |x 0 -xa| > 3 are shown in grey). The travel time distribution is ψ(s) = 1/(s α τ α + 1) 2 , α = 1, τ = 1, except in a fixed trap node α = 0.5.

Fig. 6 (

 6 Fig. 6 (Top) Example of SF network with fixed xa marked by black circle. (Bottom) FPT densities on SF network G(N, m, m 0 ) N = 100, m = 5, m 0 = 6 for different starting points x 0 and fixed xa: |x 0 -xa| = 1 (black), |x 0 -xa| = 2 (blue), |x 0 -xa| = 3 (red). The travel time distribution is ψ(s) = 1/(s α τ α + 1) 2 , α = 1, τ = 1, except in a fixed trap node α = 0.5.

Fig. 7 (

 7 Fig. 7 (Top) Example of WS network model with β = 0.2, k = 8, N = 100. (Bottom) FPT density for this network with homogeneous travel time distributions ψ(s) = 1/(sτ + 1) for different initial points x 0 , xa = 1. Different colors correspond to different shortest paths length between different x 0 and fixed xa: |x 0 -xa| = 1 (black), |x 0 -xa| = 2 (blue), |x 0 -xa| = 3 (red),|x 0 -xa| = 4 (green).

Fig. 8 (

 8 Fig. 8 (Top) The London metro map with destination Piccadilly station xa highlighted in black. (Bottom) FPT density for the London metro with travel time distributions ψ(s) = 1/(sτ + 1) except at a trap node x h at which ψx h (s) = 1/(s α τ α + 1), α = 0.5. For each starting stations x 0 the FPT density is colored according to the distance from x 0 to xa (Piccadilly station): |x 0 -xa| = 1 (black), |x 0 -xa| = 2 (blue), |x 0 -xa| = 3 (red) (FPT densities with |x 0 -xa| > 3 are shown in grey).

Fig. 9

 9 Fig. 9 The schematic representation of the short-time, intermediate-time and long-time regimes of the FPT density on various networks.

Fig. 5 (

 5 Top) Example of SF network with a fixed x a marked in black circle. (Bottom) FPT densities on SF network G(N, m, m 0 ) N = 100, m = 1, m 0 = 6 for starting points x 0 and fixed x a such that:|x 0 -x a | = 1 (black), |x 0 -x a | = 2 (blue), |x 0 -x a | = 3 (red) (FPT densities with |x 0 -x a | > 3 are shown in grey). The travel time distribution is ψ(s) = 1/(s α τ α + 1) 2 , α = 1, τ = 1,except in a fixed trap node α = 0.5. Fig.6 (Top) Example of SF network with fixed x a marked by black circle. (Bottom) FPT densities on SF network G(N, m, m 0 ) N = 100, m = 5, m 0 = 6 for different starting points x 0 and fixed x a : |x 0 -x a | = 1 (black), |x 0 -x a | = 2 (blue), |x 0 -x a | = 3 (red). The travel time distribution is ψ(s) = 1/(s α τ α +1) 2 , α = 1, τ = 1, except in a fixed trap node α = 0.5. Fig.7 (Top) Example of WS network model with β = 0.2, k = 8, N = 100. (Bottom) FPT density for this network with homogeneous travel time distributions ψ(s) = 1/(sτ +1) for different initial points x 0 , x a = 1. Different colors correspond to different shortest paths length between different x 0 and fixed x a : |x 0 -x a | = 1 (black), |x 0 -x a | = 2 (blue), |x 0 -x a | = 3 (red),|x 0 -x a | = 4 (green).

Fig. 8 (

 8 Top) The London metro map with destination Piccadilly station x a highlighted in black. (Bottom) FPT density for the London metro with travel time distributions ψ(s) = 1/(sτ + 1) except at a trap node x h at which ψx h (s) = 1/(s α τ α + 1), α = 0.5. For each starting stations x 0 the FPT density is colored according to the distance from x 0 to x a (Piccadilly station): |x 0 -x a | = 1 (black), |x 0 -x a | = 2 (blue), |x 0 -x a | = 3 (red) (FPT densities with |x 0 -x a | > 3 are shown in grey). Fig.9The schematic representation of the short-time, intermediate-time and longtime regimes of the FPT density on various networks.

in order to navigate in the network using detection of the relevant search strategies, e.g. placing heterogeneities in different parts of a network. Using the first-passage time observables one can further design network measures in order to characterize transport efficiency and spreading. The straightforward way to characterize the transport and search efficiency is to use the concept of survivability as a measure of reachability of some quantity in a network, which we plan to explore in the future. We expect the framework presented in this paper to broaden the scopes of exploration of the feature-rich networks.

Declaration

Availability of data and material

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests.

Funding

The authors acknowledge the support under Grant No. ANR-13-JSV5-0006-01 of the French National Research Agency. L.T. acknowledges support of Centre of Research and Interdisciplinarity, Paris Descartes University, France.

Authors' contributions

LT performed numerical simulations, analyzed and interpreted the analytical and numerical results of the HCTRW model on various types of networks. DG analyzed, interpreted the analytical and numerical results. Both authors contributed in writing the manuscript. Both authors read and approved the final manuscript.