Anomalous diffusion in time-fluctuating non-stationary diffusivity landscapes

Abstract

We investigate the ensemble and time averaged mean squared displacements for particle diffusion in a simple model for disordered media by assuming that the local diffusivity is both fluctuating in time and has a deterministic average growth or decay in time. In this study we compare computer simulations of the stochastic Langevin equation for this random diffusion process with analytical results. We explore the regimes of normal Brownian motion as well as anomalous diffusion in the sub- and superdiffusive regimes. We also consider effects of the inertial term on the particle motion. The investigation of the resulting diffusion is performed for unconfined and confined motion.

---

I. Introduction

The diffusion of a tracer particle is typically characterised in terms of the mean squared displacement (MSD)corresponding to the second moment of the probability density function P(x,t) to find the particle at position x at time t. When the exponent γ = 1 the law (1) describes normal Brownian diffusion, otherwise we speak of anomalous diffusion. In the latter case, the generalised diffusion coefficient K_γ has the physical dimensions cm² s^−γ, and we distinguish subdiffusion (0 < γ < 1) and superdiffusion (γ > 1) depending on the value of the anomalous diffusion exponent γ.¹

Following a surge in microscopic techniques, diffusive phenomena of passive tracer particles can now be monitored at unprecedented resolution.² Thus, for instance, the hydrodynamic backflow effects of a Brownian particle could be directly probed.³ Even more remarkable is the rapidly growing number of experimental evidence for anomalous diffusion in dense fluids⁴ as well as in living biological cells.^5–10 The motion of various endogenous and artificial tracers in live cells was shown to be subdiffusive.^11–13 However, when active dynamics such as driving by molecular motors or cytoplasmic streaming are involved, superdiffusion may also be observed.¹⁴ Massive computer simulations of pure and protein-crowded lipid bilayer membranes demonstrate transient anomalous diffusion of both lipids and proteins, the crossover to normal diffusion being delayed with increasing disorder.^15,16 In the membranes of living cells anomalous diffusion is even observed on macroscopic time scales.^17–19 As a general physical principle for anomalous diffusion various form of crowding of the environment are considered.^20–22 We note that anomalous diffusion also occurs on the level of entire organisms, such as the subdiffusion of bacteria cells in biofilms²³ or the superdiffusion of hydra or protozoa.²⁴

Several additional studies of crowded in vitro systems demonstrate the existence of non-Fickian and/or non-Gaussian motion, for instance the glassy dynamics in membrane domains,²⁵ confined diffusion of water molecules in soft environments,²⁶ polymer diffusion on nanopillar-structured surfaces,²⁷ intermittent molecular hopping on solid–liquid interfaces,²⁸ diffusion of colloidal spheres in dense crowded suspensions²⁹ and glasses,^30,31 particle diffusion in porous media with heterogeneous and position-dependent mobilities,³² and the transport of contaminants in porous and fractured geological formations.³³ Concurrently the existence of anomalous yet Brownian diffusion—a linear time dependence of the MSD (1) accompanied by significantly non-Gaussian (exponential or stretched exponential) probability density P(x,t)—was demonstrated for the motion of colloidal beads along linear phospholipid bilayer tubes,^34,35 particle dynamics in hard sphere colloidal suspensions,³⁶ and the diffusion of nanoparticles in nanopost arrays.³⁷

Brownian motion is bound to the Gaussian shape of the probability density function by the spell of the central limit theorem and thus fully characterised by the second moment (1). In contrast, anomalous diffusion dynamics is inherently non-universal, and therefore a large variety of anomalous diffusion models exists (also with non-Gaussian probability densities), depending on the exact physical circumstances defining the dynamics.^1,5,7–10 To name but a few of these anomalous diffusion processes we recall continuous time random walks with scale free trapping time distributions³⁸ and a potential additional noise source,³⁹ general trapping models,⁴⁰ correlated diffusion processes,^41,42 fractional Brownian motion⁴³ and generalised Langevin equation motion,^43–46 as well as diffusion in disordered media and on fractal structures.^47–49

Here we focus on models based on randomly and/or deterministically varying diffusion coefficients which have recently been under intense study. We show that when the diffusion coefficient varies randomly such that its distribution has a finite width, normal diffusion emerges in the long time limit. We analyse these random diffusion processes (RDPs) in terms of the MSD and time averaged MSD (typically evaluated in single particle tracking and simulations studies), for both unconfined and confined diffusion. In addition, we quantify the degree of randomness between different individual trajectories.

In Section II we provide a concise overview of heterogeneous diffusion processes, followed by a definition of the various observables in Section III. In Section IV we present the details of the specific model investigated here and the numerical scheme used to simulate the RDPs. We then present the main results of our calculations for the ensemble averaged MSD, time averaged MSD, probability distribution function P(x,t), and the ergodicity breaking parameter of RDPs in Sections V and VI, respectively, for massless and massive particles. We consider the situation both in the absence and in the presence of an external confinement. Section VII summaries our findings and discusses their possible applications and generalisations.

II. Heterogeneous mobility: models and examples

A. Random diffusivity models

For massless particles the study of normal and anomalous diffusion in the presence of random diffusivity fields recently attracted considerable attention.^19,50–53 Several models assume that the instantaneous diffusion coefficient D is governed by a steady state distribution p(D) in an annealed fashion, that is, the instantaneous value of D is independent of the actual particle position x. Particular attention received the idea of a diffusing diffusivity introduced by Chubynsky and Slater.⁵¹ They assume an exponential distribution

p(D) = D₀⁻¹ exp[−D/D₀] (2)

and weigh the standard Gaussian P(x,D,t) = (4πDt)^−1/2 exp(−x²/[4Dt]) with this function, to obtain the exponential probability distribution function (PDF)The characteristic decay length of this PDF grows with the diffusion time as ≃t^1/2,⁵¹ but the MSD grows linearly with time, that is, follows eqn (1) with γ = 1 and K_γ = D₀.⁵¹ Other than exponential diffusivity distributions—for instance, power-law forms of p(D)—were shown to lead to subdiffusive and non-ergodic MSD behaviour.⁵⁰ After the current manuscript was submitted, the authors became aware of the simulation-based study⁵⁴ of particle diffusion in rough energy landscapes with both Gaussian and Gamma distributed local energy values.

Dynamical heterogeneities—as reflected in the above assumption of a random diffusivity—are considered a characteristic property of systems such as supercooled or glassy liquids.^52,55–58 Quite broad distributions of particle diffusivities were detected in a number of living systems, for instance, for the motion of pathogen receptors on two dimensional cell membranes¹⁹ (see also ref. 59), the motion of Cajal bodies in eukaryotic nuclei,⁶⁰ one dimensional diffusion of repressor proteins on the DNA,⁶¹ and for the motion of proteins along the corrugated landscape created by the DNA sequence.^62,63 In these systems, the inherent stochasticity of the diffusive properties of a tracer particle as well as the heterogeneities of its environment contribute to the observed distribution of diffusivities subsumed in the distribution p(D).

The reader is particularly referred to the characterisation of the dynamical spreading of a population of nematode worms, both in homogeneous and heterogeneous environments.⁶⁴ Other examples of living systems with diffusing individuals obeying non-Gaussian distributions of diffusivities, speeds of motion v, or turning angles are also mentioned in this study. Hapca et al.⁶⁴ state that the anomalous diffusion monitored in the heterogeneous populations of worms can be solely due to fat-tailed, e.g., Gamma distributed⁶⁵ forms of diffusivities, while the motion of each individual remains Brownian. This study served as a strong biological motivation for us in trying to unveil the properties of particle diffusion with a given time dependent form of p(D,t).

Also note that fat-tailed leptokurtic distribution of particle mobilities—often occurring in population of individuals of a species—can ensure a facilitation of their colonial invasion,⁶⁴ as compared to the standard Brownian diffusion law of spreading. Such skewed distributions p(D) or p(v) (of the particle speeds) can originate from medium heterogeneities when the organisms explore different regions of space with different mobilities. The notion of fat-tailed distributions and faster than standard front propagation emerges also in long-distance dispersal of plants^66,67 and pollen,⁶⁸ in patterns of fish movements,⁶⁹ as well as in rare event driven spreading of plant pathogens.⁷⁰

B. Deterministic variation of the diffusivity with position or time

Following experimental observations of deterministic gradients of the local diffusivity in both pro- and eukaryotic cells^71,72 and the existence of thermal gradient conditions,⁷³ the model of heterogeneous diffusion processes (HDPs) with a power-law, exponential, and logarithmic form for D(x) was recently introduced by the authors;^74–76 see also ref. 77–80. These Markovian processes based on a Langevin description with multiplicative noise exhibit anomalous diffusion and weak ergodicity breaking.^74–76 The latter emerges due to the fact that even in the limit of long trajectories time and ensemble averages of physical observables do not coincide,^6,8 see below. Models with a power-law time dependence

D(t) ≃ t^α−1 (4)

of the diffusivity exist, the so-called scaled Brownian motion (SBM).^8,81–87 SBM was originally introduced by Batchelor in the description of Richardson turbulence.⁸⁸ Note that in the context of such highly non-stationary processes the degree of ergodicity breaking is controlled via introducing a time- or length-dependent scale into the problem.^45,89 A combined space-time diffusivity dependence of the form D(x,t) ≃ |x|^βt^α−1 was also investigated.^90–92 Aged and confined versions of these processes were recently considered as well.^76,93 Here ageing refers to the explicit dependence of the process on the overall time span of its evolution, as expected for non-stationary processes. The limiting cases with scaling exponent β → 2 for HDPs and α → 0 for SBM were shown to lead, respectively, to an ultrafast (exponential) and ultraslow (logarithmic) MSD growth with time.^94–96

C. Massless versus massive particles

For massive particles, the situation in anomalous diffusion is often less clear. While the underdamped limit of standard Brownian motion⁹⁷ and of fractional Langevin equation motion⁴⁶ are well understood, for other anomalous diffusion processes these limits and general solutions are just emerging. In particular, for SBM it was recently shown that the long time limit of underdamped motion (including the inertia term) in general does not correspond to the overdamped limit of the same motion.⁹⁸ Also, the recent study⁹⁹ addresses a giant particle diffusion in the underdamped limit with a temperature dependent diffusion coefficient and in the presence of a bias.

One purpose of the current study is to investigate RDPs for massive and massless particles. A growing interest in the diffusive behaviour of tracked particles combined with the unprecedented precision of experimental observations, in particular, at short times for the diffusion of small particles in living cells,¹⁰⁰ pose a need for the development of new and more flexible models of stochastic processes. Thus, a larger pool of theoretical models is necessary for quantitative descriptions of these systems, with possibly fewer number of model parameters.

Some implications of a finite particle mass for diffusion processes with position dependent diffusivity of the form D(x) ≃ |x|^β were recently examined.¹⁰¹ The reader is also referred to the studies^102,103 regarding the inertial Langevin dynamics in media with space inhomogeneous friction, and conventions of how to interpret the associated multiplicative stochastic equation as well as the existence of fluctuation–dissipation relations for such systems. In what follows, we refer to the diffusion coefficient D as to local variable in space and time, rather as to a long time asymptote of the Einstein relation, see the discussion in ref. 103.

III. Observables of diffusion processes

Anomalous diffusion processes can be classified by the MSD diffusion exponent γ. In single particle tracking and simulations studies garnering few but long individual time series x(t) of the particle position the time averaged MSD^6,8is typically employed. Here T is the total length of the trajectory (observation time) and Δ is the lag time. Note that while the ensemble averaged MSD (1) is a spatial average at a particular time instant t, the time averaged MSD (5) for any given lag time Δ is taken over the entire history of the trajectory x(t). As usual, ensemble averaging is denoted hereafter by angular brackets, while time averaging is indicated by the overline. To obtain smoother curves for the time averaged MSD an additional average is taken over N trajectories, defining the mean time averaged MSD^6,8

Ergodicity in the Boltzmann–Khinchin sense typically assumed in equilibrium statistical mechanics would imply the equivalence of ensemble and time averaged MSD in the limit of long measurement times, . Following Bouchaud¹⁰⁴ the breakdown of this relation is referred to as weak ergodicity breaking,^{5,8,105–108}

Continuous time random walks and HDPs are known to be weakly non-ergodic,⁸ while diffusion on fractals is ergodic on the infinite cluster but not on the entirety of all clusters.⁴⁹ In contrast, other diffusive processes such as fractional Brownian motion and SBM are only marginally non-ergodic.^{8,81,86,87,109–113}

A distinctive measure of non-reproducibility of individual time averaged MSD traces is the ergodicity breaking parameter¹⁰⁸

EB(Δ) = 〈ξ²(Δ)〉 − 1 (8)

based on the dimensionless ratio quantifying the spread of individual time averaged MSDs about their mean (6). Typically EB of a weakly non-ergodic process decays to zero with increasing trace length slower than for the standard Brownian motion,^46,112Or, EB may even attain a finite value as Δ/T → 0, for instance, for HDPs and continuous times random walks.^{8,74,106,108,114} Often, also the ratio of the time and ensemble averaged MSDs¹¹⁵provides additional information about the ergodic properties of the diffusion process.

IV. The random diffusivity model

In this section we describe the details of RDPs. As a generalisation of eqn (2), the instantaneous value of the diffusion coefficient on each simulation step is independently chosen from the Rayleigh distributionThe mean particle diffusivity is then given byand the diffusivity variance is 〈(D− 〈D〉)²〉 = (2 − π/2)D_σ². Independence of successive values of the diffusion coefficient indicates no temporal correlations in its fluctuations. Our system is thus out of equilibrium (the temperature is not fixed) and the fluctuation dissipation theorem does not hold.

The PDF (11) is a smooth function in the range from D = 0 to D = ∞ and it vanishes on the boundaries of this interval. This distribution p(D) is used instead of a Gaussian distributed diffusion coefficient to avoid non-physical negative D values. The distribution (11) is thus physically different from the exponential p(D) form given by eqn (2) used by Chubynsky and Slater.⁵¹ When the mean diffusivity stays constant over time, in the long time limit the particle diffuses normally. However, when the mean diffusivity is allowed to vary in addition deterministically as a powerlaw,

〈D(t)〉 ≃ t^ω, (13)

the resulting process is reminiscent of SBM. Physically, such an increase of the mean diffusivity could be due to a diffusing diffusivity of the form 〈D(t_i)〉 = |〈D(t_i−1)〉 + ζ_D(t_i−1)|, where ζ_D(t_i) is an incremental change. This is analogous to the power-law growth of the waiting times in the correlated continuous time random walks.⁴²

We consider below both massless and massive particles diffusing in both unconfined and a confined environments. We implement the same algorithms for the iterative computation of the particle displacement x(t) as developed for HDPs⁷⁴ and combined HDP-SBM motion.⁹² First, we simulate the one dimensional overdamped Langevin equation

driven by zero-mean and unit-variance Gaussian noise ζ(t). At step i + 1 the particle displacement is given bywhere (y_i+1 − y_i) are the increments of the Wiener process. Unit time intervals separate consecutive steps. To avoid possible particle stalling we regularise D by adding a small constant D₀ = 10⁻³.^75,76 This does not affect the intermediate- and long-time diffusive behaviour. The particle's initial position is x₀ = x(t = 0) = 0.1. In the second part of the paper, we simulate the underdamped Langevin equation for a particle of mass m,with the unit damping coefficient set below to η = 1. At a step i + 1 the particle displacement is found from the iteration schemewhere the instantaneous diffusivity is taken from eqn (11) with the mean (13).

V. Results: overdamped motion

A. Free diffusion

We start with the diffusion of massless particles with a fluctuating diffusivity and time invariant mean. As naively expected, we find that due to friction the MSD in the long time limit is Brownian,

〈x²(t)〉 ∼ x₀² + 2〈D(t)〉t. (18)

The MSD and the time averaged MSD are nearly identical after a fast relaxation of the starting position x₀ of the particle and the ergodicity is approximately fulfilled at all times (results not shown here).

Now we address the more interesting case of RDPs with instantaneous diffusivity chosen from the distribution (11) with a time dependent mean (13). Namely, the most likely diffusivity D_σ(t) at simulation step i is

where the prefactor is chosen for convenience and the coefficient [scr D, script letter D] tunes the magnitude of the diffusivity. From eqn (14) we straightforwardly compute

〈x²(t)〉 = x₀² + 2[scr D, script letter D]t^γ, (20)

with

γ = 1 +ω. (21)

Fig. 1 depicts the case of ω = 1/2 or the MSD diffusion exponent γ = 3/2. As can be seen from the simulations the time averaged MSD grows linearly with the lag time Δ, as in the Brownian case. This process also reveals a quite moderate amplitude spread of individual traces , see the thin red curves in Fig. 1. Obviously when the lag time approaches the observation time, Δ ∼ T, the amplitude scatter increases due to the deteriorating statistic of .⁸ Varying [scr D, script letter D] in Fig. 1 we demonstrate that, as expected, larger initial diffusivities give rise to a faster approach of the MSD to the theoretical asymptote. In contrast, for rather small diffusivities (smaller [scr D, script letter D] values) the system needs more time to approach the long time asymptote, Fig. 1. A diminished magnitude of the MSD at smaller [scr D, script letter D] values inevitably leads to a decrease in the magnitude of the time averaged MSD, see below. At intermediate to long times the superdiffusive MSD regime with γ = 3/2 emerges. Finally, towards the very end of the trace the MSD and the time averaged MSD coincide, as they should.⁸

Fig. 1 MSD 〈x²(t)〉 and time averaged MSD (thick blue curves) as well as individual time traces (red curves) for overdamped RDPs. The asymptotes (20) and (22) for the MSD and the time averaged MSD are shown by the dashed curves. The asymptotes often superimpose with the results of simulations. Parameters: the scaling exponent is ω = 1/2, the trace length is T = 10⁵, and number of traces used for the averaging is N = 150. The starting position is x₀ = 0.1 and the parameter takes the values 10⁻¹, 10⁻², 10⁻³, and 10⁻⁴ for the panels from left to right.

Analytically, we obtain for the time averaged MSD from eqn (14) for x(t) and after averaging over the noise ζ(t) and diffusion coefficient realisations p(D) the result

This expression nicely agrees with the results of computer simulations for all ω and [scr D, script letter D] values investigated, see Fig. 1 and 8. It is not surprising that both MSD and time averaged MSD are proportional to [scr D, script letter D] determining the basal value of the particle diffusivity, eqn (19). In the limit of short lag times, Δ ≪ T, from eqn (22) we recover the scaling behaviourThus, for superdiffusive RDPs with γ > 1 the magnitude of the time averaged MSD is a growing function of the trace length T, while for subdiffusive RDPs magnitude decreases with T, in agreement with Fig. 9. In the limit Δ→T the ensemble and time averaged MSDs coincide,⁸ as it is easy to check from eqn (20) and (22) and corroborated in Fig. 1. The non-equivalence of the time averaged MSD (22) and the time averaged MSD (20) demonstrates that the system is weakly non-ergodic. The scaling behaviour (23)—regarding the magnitude of the time averaged MSD in terms of the power law of the trace length T and the linearity in the lag time—is analogous to that obtained for subdiffusive continuous time random walks^{13,106,108,116} and their correlated version⁴² as well as for HDPs^74,75,92 and SBM;^87,93 see also ref. 8 for an overview.

For RDPs the ergodicity breaking parameter EB computed from simulations tends to follow the asymptote (9) for Brownian motion at intermediate and long times, Fig. 2. We observe that the initial relaxation of EB to this asymptote is relatively fast. As we show in Fig. 2, after this relaxation time the parameter (10) becomes a power-law function of the lag time Δ,

Fig. 2 Ergodicity breaking parameter EB and the parameter versus lag time Δ for overdamped RDPs. The curves are computed for the parameters of Fig. 1 with ω = 1/2. The curves for from top to bottom correspond to , respectively. The asymptote (9) of Brownian motion and the relation (24) are the dashed lines.

Diffusivity distributions p(D) whose mean diffusivity grows with time may be viewed to correspond to an effectively increasing temperature in the system. The opposite case of a temporally shrinking width may stem from a cooling of the system in the course of time. In this respect the current process is reminiscent of SBM.^81,87,93 An important example for the latter are granular gases with a relative velocity dependent restitution coefficient.⁸⁵

The particle spreading for very subdiffusive RDPs can be compared to the PDF of SBM, identical to that of fractional Brownian motion for x₀ = 0 for natural boundary conditions of a vanishing PDF of diffusing particles at |x| → ∞.^83,87 Namely, after the substitution of the corresponding MSD (20) this produces

Fig. 10 compares the conjecture (25) with the result from simulations of RDPs. We see, however, that as the MSD scaling exponent increases a distinct cusp of the PDF starts to develop at the origin. The PDF of superdiffusive RDPs becomes pronouncedly non-Gaussian. This spike however cannot be described by a convolution of the diffusivity distribution (11) with the kernels of Brownian motion and SBM motion. A more detailed investigation is required to understand this spike at short times and possibly exponential forms of the PDF tails at long times. This generally non-Gaussian and γ-dependent shape of the PDF is one important distinction of RDPs with time dependent and in addition fluctuating diffusivities, as compared to the SBM process with the deterministic value of the diffusion coefficient at each step, D(t) = 〈D(t)〉.

B. Confined diffusion

We now turn to confined RDPs on an interval −L < x < L. Such confined motion is important especially for the understanding of diffusion processes in biological cells. In cells—due to their external confinement by the plasma membrane and internal compartmentalisation—a diffusing tracer frequently collides with boundaries. As expected, after an initial free diffusion the MSD converges to the stationary plateau⁸as demonstrated in Fig. 3. The time averaged MSD, by virtue of its definition (5), approaches twice the value of the MSD in the long time limit.

The existence of a plateau is similar to that of standard Brownian motion, and interval-confined SBM and HDPs.^76,92 Note that SBM confined by an external potential has a time dependent thermal value of the MSD.^87,93 The behaviour of confined continuous time random walks is strikingly different, there confinement leads to a crossover to a second power-law regime in the time averaged MSD.^13,116 The PDF of confined RDPs approaches a uniform distribution of particles on the interval.

Fig. 3 Same as in Fig. 1 but for confined overdamped RDPs. The box size is L = 10², ω = 1/2, and . The asymptote (26) and twice this value are the dashed lines.

At more severe confinement the ergodicity breaking parameter EB at large lag times Δ values starts to deviate form the Brownian asymptote (9), see Fig. 4. In addition, we find that for a fixed width of the confining interval and varying trace length T the EB parameter follows the scaling relation

EB(T) ≃ 1/T, (27)

as illustrated in Fig. 11A and B. This is a standard decrease of the EB parameter for longer trajectories, a property ubiquitous among a number of both ergodic and non-ergodic stochastic processes.⁸ The decrease of EB with T indicates a progressively more ergodic diffusion for longer particle traces.

Fig. 4 Parameters EB and for confined overdamped RDPs for the parameters of Fig. 3. The Brownian asymptote (9) and the plateau are shown as the dashed lines.

VI. Results: underdamped motion

A. Free diffusion

In this section, we study the diffusion of massive particles in the same time dependent random diffusivity scenario (19) based on the underdamped Langevin equation. In particular, we explore to what extent inertia effects modify the long time behaviour of the MSD and the time averaged MSD, as compared to the overdamped RDPs considered above. The general solution for the particle MSD follows from the standard procedure for the Brownian motion of massive particles.⁹⁷ Namely, we obtainwhere we defined θ = η/m and A² = 2γ[scr D, script letter D]/m². Moreover, v₀ is the initial particle velocity. Then the MSD of the particles after averaging over the noise ζ can formally be written asHere denotes the generalised exponential integral.

For zero initial velocity of the particles v₀ = 0, as in the computer simulations performed here, the inertial term in the Langevin equation gives rise to the initial MSD scaling of the form 〈x²(t)〉 ∼ t^γ+2. It is due to progressively accelerating (heating) particles. Explicitly, for the MSD at short times we get

This faster than ballistic MSD regime often called hyperdiffusion is known to emerge, for instance, for a power-law like transient heating of particles with temperature variation of the form .¹¹⁷ This superballistic behaviour emerges for RDPs with fluctuating diffusivities and time dependent mean, in analogy with a faster than linear short time ballistic regime in Brownian motion.⁹⁷

Note that the short time MSD regime (30)—with the scaling exponent by one larger than the long time MSD exponent—is absent in the model of underdamped SBM elucidated by us recently in ref. 85. The reason is that the damping coefficient is set to be temperature independent in the current model, whereas in the model of underdamped SBM η(t) is coupled to the diffusivity variations via the generalised time-local Einstein relation⁸⁵

D(t) = k_BT(t)/(mη(t)). (31)

So, the fluctuation–dissipation theorem is valid, contrary to the current approach. For the underdamped SBM process, the relation η(t) ∼ T^1/2(t) is consistent with the physical picture of elastically colliding and relaxing particles in a bath with a deterministically varying temperature.⁹⁸ The reader is also referred to ref. 118 for studying different relationships between the friction coefficient and velocity for passive and active¹¹⁹ particles, including nonlinear forms.

Note that the diffusive and ergodic properties of underdamped SBM were recently considered as well.⁹⁸ It was demonstrated that inertial effects relax rather quickly in the course of particle diffusion for superdiffusive (γ > 1) situations. This follows from comparing the magnitudes of the acceleration and friction terms in the Langevin equation. Conversely, for small positive values of γ a finite particle mass yields an extensive intermediate regime, both for the MSD and the time averaged MSD growth behaviour with time. Interestingly, in the case of ultraslow logarithmic SBM motion—realised for the limiting value γ = 0—the overdamped limit of particle diffusion⁹⁶ is not reached at long times, independent on the total measurement time.⁹⁸

The time averaged MSD of underdamped RDPs follows from eqn (28) and (5),

This integral expression can be evaluated numerically. In the limit of short lag times Δ ≪ T we can evaluate the integral and findHere is the generalised incomplete Gamma function and Γ(a,0) = Γ(a) is the Gamma function. The short time asymptotes of both the MSD and the time averaged MSD are plotted as the dot-dashed curves in Fig. 5 showing nice agreement with the results of computer simulations of eqn (15).

Fig. 5 Ensemble averaged MSD, individual time averaged MSDs, and mean time averaged MSD, with the same notation as in Fig. 1, evaluated for underdamped RDPs with particle mass m = 1, 5, 10, and 50, for panels from left to right. N ∼ 10² traces are shown. The long time asymptotes (20) and (22) are the dashed lines. The dot-dashed lines are the asymptotes for the short time regimes, eqn (30) and (34). Parameters: , ω = 1/2, and T = 10⁵.

Expanding the Gamma functions in the corresponding limit we find from expression (33) for light particles or high friction in the system—that is for θT ≪ 1—that

This has the form of the short time MSD behaviour of standard Brownian motion.⁹⁷ All particle masses in this study guarantee the validity of eqn (34) as the short time expansion for the time averaged MSD in this underdamped limit. Comparing the -variation with eqn (34) in Fig. 9B further supports this validity regarding the dependence on the length of particle trajectory T. In the opposite limit of θT ≪ 1—very massive particles or low friction for the particle motion—we find

Let us now describe the results of our computer simulations and compare them with these analytical predictions. We observe that in the particle displacement the initial condition of zero particle velocity (v₀ = 0) relaxes within several initial diffusion steps. Naturally, it takes for the system longer to accelerate heavier particles, as shown in Fig. 5, in agreement with eqn (30). The initial slow acceleration of heavy particles yields slowly growing MSD that gets in turn reflected in small amplitudes of the time averaged MSD at short lag times Δ. The MSD follows eqn (30) for short times and then crosses over to the long time scaling (20). Note that the same initial quadratic regime was observed in ref. 98 for the short time behaviour of the time averaged MSD of the standard underdamped SBM process. For the time averaged MSD the heavier particles feature a longer ballistic regime, see Fig. 5. The magnitude of decreases with the particle mass, in agreement with eqn (34). At longer times the MSD and time averaged MSD approach the results expected for the overdamped RDP motion, shown as the long time asymptotes in Fig. 5. As expected, apart from the initial ballistic regime of the time averaged MSD described by eqn (34), the analytical solution for the overdamped limit (22) describes the long time behaviour of the results of our computer simulations.

Similar to the overdamped situation, at short lag times the magnitudes of follow the relation (23), see Fig. 9B. We also observe that for underdamped RDPs the spread of individual traces is similar to that of standard Brownian motion. This small spread is consistent with the observation that the EB parameter for this underdamped RDP motion does not deviate strongly from the Brownian asymptote (9) at intermediate and long times, see Fig. 6. In the region of short Δ the deviations are quite substantial, particularly for massive particles exhibiting a ballistic initial growth of and a nonlinear growth of the MSD, see Fig. 5. The auxiliary parameter —again after the initial particle acceleration—follows the asymptote (24), see the thick curves in Fig. 6. The deviations of EB and at short lag times from the Brownian asymptote is more evident for massive particles.

Fig. 6 Parameters EB and versus lag time Δ for underdamped RDPs. The Brownian asymptote (9) and eqn (24) are the dashed lines. The curves correspond to varying particle mass, m = 1, 5, 10, 50; the other parameters are the same as in Fig. 5.

B. Confined diffusion

We complete the analysis of RDPs with the study of the underdamped motion in a confining box. We observe that at short time the MSD develops similar to the unconfined scenario. Once the boundary of the confined region is reached, the plateaus start to develop at the same levels as for the overdamped RDP case both for the MSD and time averaged MSD, see eqn (26) and Fig. 7. Towards the very end of the trajectory at Δ = T, the MSD and the time averaged MSD coincide, as they should.⁸ The longer the entire trajectory, however, the narrower the range of lag times where this convergence takes place, and thus the more precise should be the Δ-sampling in this region—that is often computationally costly. This effect was studied in detail for the pure SBM motion confined in harmonic potentials⁸⁷ and for HDPs confined between hard walls.^76,92

Fig. 7 Analogue of Fig. 5 for confined underdamped RDPs with particle mass m = 1, 5, 10, and 50, from left to right. We show N = 150 traces for each set of parameters and the interval width is L = 10². The asymptotes for the confined motion (26) are the dashed lines.

Fig. 8 Same as in Fig. 1 but at ω = −1/2 (A) and ω = 1 (B) computed for and N = 150. The asymptotes for the MSD and time averaged MSD are shown as the dashed lines.

Fig. 9 Time averaged MSD versus trace length T for (A) non-confined overdamped RDPs for ω = −1/2, 1/2, 1 (for data points from bottom to top) and for (B) non-confined underdamped RDPs for m = 1, 5, 10, 50 (for data points from top to bottom) and ω = 1/2. Other parameters are the same as in Fig. 1 and . Dashed lines indicate the scaling relation (23) in panel (A) and eqn (34) in panel (B).

Fig. 10 PDF of the overdamped RDPs for ω = −1/2, 1/2, and 1 (from left to right). Other parameters are the same as in Fig. 1 with and N = 300. The PDF asymptote for subdiffusive SBM (25) is shown as the dashed curve.

Fig. 11 EB parameter at Δ = 1 computed for free and confined overdamped RDPs (panels (A) and (B)) and for free and confined underdamped RDPs (panels (C) and (D)). Parameters: for all the panels ω = 1/2, and (panel (A): free overdamped process) m = 0, L = ∞; (panel (B): confined overdamped process) m = 0, L = 10²; (panel (C): free underdamped process) m = 50, L = ∞; (panel (D): confined underdamped process) m = 50, L = 10².

The PDFs of confined underdamped RDPs at varying box width L is presented in Fig. 12. We observe that for wide intervals the particles are nearly uniformly distributed on the interval (see the dashed lines in Fig. 12), with only insignificant increase in the particle occupancies near the box boundaries due to reflections. Note that the particle starting position at x = x₀ is still slightly visible in the PDF for a weak confinement. As the confinement becomes more severe, the particle accumulation near the interval boundaries occupies a larger fraction of space available for diffusion.

Fig. 12 PDF of underdamped confined RDPs for L = 3, 10, 30, and 100 (from top to bottom), computed for ω = 1/2, , T = 10⁵, m = 50. Dashed lines designate the uniform distribution of particles on the interval.

The EB parameter for massive particles deviates progressively from the Brownian law (9) at both short and long lag times Δ (not shown). This is due to slow particle acceleration at short times (MSD plateau) and particle confinement at long times, respectively. For free and confined underdamped RDPs—similarly to the overdamped situation—the EB parameter follows the asymptote (27) with the trajectory length T, see Fig. 11C and D. As the confinement becomes less severe, the EB parameter approaches the value for the free underdamped RDP motion. Fig. 13 illustrates this EB evolution with the width of the confining interval.

Fig. 13 EB parameter at Δ = 1 for underdamped RDPs with varying degree of confinement L. Parameters are the same as in Fig. 12: ω = 1/2, , T = 10⁵, m = 50.

VII. Conclusions

We examined the ensemble and time averaged characteristics of random diffusion processes. The randomness of the diffusion coefficient D was implemented in the model via a non-stationary distribution p(D). RDPs are not thermalised, that is, the motion of the particles is inherently out of equilibrium. The distribution of the diffusion coefficient reflects individual variations of particle diffusivities and heterogeneities of the environment. For typical out-of-equilibrium systems such as biological cells this does not pose any restrictions to our model. RDPs represent a quite flexible model to study asymptotically Brownian and anomalous diffusive systems with a locally fluctuating diffusion coefficient to model physical situations in many complex systems.

We rationalised both by computer simulations and analytically the MSD, the time averaged MSD, and the ergodicity breaking parameter of RDPs. Unconfined and interval-confined motion was examined. We found that in terms of these standard characteristics subdiffusive RDPs appear similar to subdiffusive SBM with a deterministic diffusivity variation in time. For superdiffusive RDPs the fluctuations grow with time. Concurrently, the average diffusion coefficient 〈D(t)〉 grows with time together with the spread of its values. These features reflect an increasing temperature and more pronounced fluctuations of the medium in the course of particle diffusion. The properties of ageing overdamped and underdamped RDPs will be considered elsewhere and compared to the ageing properties of SBM,⁹³ HDPs,⁷⁶ and continuous time random walks.¹¹⁴

Living cells feature heterogeneous and densely crowded environments established by a “melange” of various macromolecules and (importantly) a rather viscoelastic solution between them. This often leads to a broadening in the distribution of diffusion coefficients and subdiffusive exponents, as observed for obstructed diffusion of various tracers.^{12,59,120,121} In particular, some extensions of the standard diffusion models to account for these effects—similar to our p(D) distribution for the SBM like diffusion model presented above—appear necessary e.g. for a quantitative fit of fluorescence recovery after photobleaching curves.¹²⁰ The models of SBM type—with the diffusivity formally decaying in time according to the power law (13)—are often implemented to describe the subdiffusive MSD behaviour (1) of the tracer particles in cells. This anomalous MSD scaling was observed e.g. via fitting the shape of the autocorrelation curves of fluorescence correlation spectroscopy measurements.^12,122,123

In single particle tracking measurements in biological cells some time dependent scatter of the diffusion coefficients was also detected.^59,124 It is necessary in theoretical models i.a. to distinguish between the normal, restricted, and fully trapped populations of the tracers. In fact, a Gamma distribution similar to the Rayleigh distribution (11) used above was proposed in ref. 59 to characterise the scatter of the MSD distribution of the diffusing particles. On the level of diffusing simple organisms, Gamma like diffusivity distributions were documented for the motion of nematode worms.⁶⁴ The latter also exhibit non-Gaussian PDFs of the particle displacements with a “spike” at the origin,^64,125 similar to some of our findings. Also, the recent study¹²⁶ of anomalous and non-ergodic dynamics of particles within a predator-prey model with a broad distribution of diffusion coefficients of interacting partners needs to be mentioned here.

The current study with its preset functional form of the diffusivity distribution and a deterministic law (13) represents a first step into the terrain of stochastic processes with fluctuating and time varying diffusivities. A more general consideration would correspond to a system of coupled stochastic differential equations for the particle position and its diffusivity. The first equation is the standard Langevin equation, while the second equation involves an additional, generally decoupled noise source governing D(t) variation. The correlation function and other noise properties—not necessarily Gaussian—determine then both the ensemble and the time averaged MSDs of diffusing particles.¹²⁷

Appendix

In this Appendix we present several additional figures (Fig. 8–13) supporting the claims in the main text of the manuscript.

References

1. J.-P. Bouchaud and A. Georges, Phys. Rep., 1990, 195, 127; R. Metzler and J. Klafter, Phys. Rep., 2000, 339, 1; R. Metzler and J. Klafter, J. Phys. A: Math. Gen., 2004, 37, R161.
2. C. Bräuchle, D. C. Lamb and J. Michaelis, Single Particle Tracking and Single Molecule Energy Transfer, Wiley-VCH, Weinheim, Germany, 2012; X. S. Xie, P. J. Choi, G.-W. Li, N. K. Lee and G. Lia, Annu. Rev. Biophys., 2008, 37, 417.
3. T. Franosch, M. Grimm, M. Belushkin, F. M. Mor, G. Foffi, L. Forro and S. Jeney, Nature, 2011, 478, 7367.
4. D. S. Banks and C. Fradin, Biophys. J., 2005, 89, 2960; G. Guigas, C. Kalla and M. Weiss, Biophys. J., 2007, 93, 316; J. Szymanski and M. Weiss, Phys. Rev. Lett., 2009, 103, 038102; W. Pan, L. Filobelo, N. D. Q. Pham, O. Galkin, V. V. Uzunova and P. G. Vekilov, Phys. Rev. Lett., 2009, 102, 058101; J.-H. Jeon, N. Leijnse, L. B. Oddershede and R. Metzler, New J. Phys., 2013, 15, 045011.
5. S. Burov, J.-H. Jeon, R. Metzler and E. Barkai, Phys. Chem. Chem. Phys., 2011, 13, 1800.
6. E. Barkai, Y. Garini and R. Metzler, Phys. Today, 2012, 65(8), 29.
7. F. Höfling and T. Franosch, Rep. Prog. Phys., 2013, 76, 046602.
8. R. Metzler, J.-H. Jeon, A. G. Cherstvy and E. Barkai, Phys. Chem. Chem. Phys., 2014, 16, 24128.
9. Y. Meroz and I. M. Sokolov, Phys. Rep., 2015, 573, 1.
10. C. Manzo and M. F. Garcia-Parajo, Rep. Prog. Phys., 2015, 78, 124601.
11. A. Caspi, R. Granek and M. Elbaum, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 011916; G. Seisenberger, et al. , Science, 2001, 294, 1929; I. Golding and E. C. Cox, Phys. Rev. Lett., 2006, 96, 098102; I. Bronstein, Y. Israel, E. Kepten, S. Mai, Y. Shav-Tal, E. Barkai and Y. Garini, Phys. Rev. Lett., 2009, 103, 018102; S. M. A. Tabei, S. Burov, H. Y. Kim, A. Kuznetsov, T. Huynh, J. Jureller, L. H. Philipson, A. R. Dinner and N. F. Scherer, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 4911.
12. M. Weiss, M. Elsner, F. Kartberg and T. Nilsson, Biophys. J., 2004, 87, 3518.
13. J.-H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber-Unkel, K. Berg-Sorensen, L. Oddershede and R. Metzler, Phys. Rev. Lett., 2011, 106, 048103.
14. A. Caspi, R. Granek and M. Elbaum, Phys. Rev. Lett., 2000, 85, 5655; N. Gal and D. Weihs, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 020903(R); D. Robert, T. H. Nguyen, F. Gallet and C. Wilhelm, PLoS One, 2010, 4, e10046; J. F. Reverey, J.-H. Jeon, H. Bao, M. Leippe, R. Metzler and C. Selhuber-Unkel, Sci. Rep., 2015, 5, 11690.
15. G. R. Kneller, K. Baczynski and M. Pasienkewicz-Gierula, J. Chem. Phys., 2011, 135, 141105; T. Akimoto, E. Yamamoto, K. Yasuoka, Y. Hirano and M. Yasui, Phys. Rev. Lett., 2011, 107, 178103; J.-H. Jeon, H. Martinez-Seara Monne, M. Javanainen and R. Metzler, Phys. Rev. Lett., 2012, 109, 188103; M. Javanainen, H. Hammaren, L. Monticelli, J.-H. Jeon, R. Metzler and I. Vattulainen, Faraday Discuss., 2013, 161, 397; S. Stachura and G. R. Kneller, J. Chem. Phys., 2014, 40, 245; E. Yamamoto, A. C. Kalli, T. Akimoto, K. Yasuoka and M. S. P. Sansom, Sci. Rep., 2015, 5, 18245.
16. J.-H. Jeon, M. Javanainen, H. Martinez-Seara, R. Metzler and I. Vattulainen, Phys. Rev. X, 2016, 6, 021006.
17. A. V. Weigel, B. Simon, M. M. Tamkun and D. Krapf, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 6438; D. Krapf, Curr. Top. Membr., 2015, 75, 167; D. Krapf, G. Campagnola, K. Nepal and O. B. Peersen, Phys. Chem. Chem. Phys., 2016, 18, 12633.
18. R. Metzler, J.-H. Jeon and A. G. Cherstvy, Biochim. Biophys. Acta, 2016, 1858, 2451.
19. C. Manzo, J. A. Torreno-Pina, P. Massignan, G. J. Lapeyre, Jr., M. Lewenstein and M. F. Garcia Parajo, Phys. Rev. X, 2015, 5, 011021.
20. H. Berry and H. A. Soula, Front. Physiol., 2014, 5, 437; H. Berry and H. Chate, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 022708; F. Trovato and V. Tozzini, Biophys. J., 2014, 107, 2579; M. Weiss, Int. Rev. Cell Mol. Biol., 2014, 307, 383; M. J. Saxton, J. Phys. Chem. B, 2014, 118, 12805; D. S. Banks, C. Tressler, R. D. Peters, F. Höfling and C. Fradin, Soft Matter, 2016, 12, 4190; T. Sentjabrskaja, E. Zaccarelli, C. De Michele, F. Sciortino, P. Tartaglia, F. Höfling and T. Franosch, Phys. Rev. Lett., 2007, 98, 140601; S. Leitmann and T. Franosch, Phys. Rev. Lett., 2013, 111, 190603; T. Voigtmann, S. U. Egelhaaf and M. Laurati, Nat. Commun., 2016, 7, 11133.
21. S. K. Ghosh, A. G. Cherstvy and R. Metzler, Phys. Chem. Chem. Phys., 2015, 17, 472.
22. S. Ghosh, A. G. Cherstvy, D. S. Grebenkov and R. Metzler, New J. Phys., 2016, 16, 013027.
23. S. S. Rogers, C. van der Walle and T. A. Weigh, Langmuir, 2008, 24, 13549.
24. A. Upadhyaya, J.-P. Rieu, J. A. Glazier and Y. Sawada, Physica A, 2001, 293, 549; L. G. A. Alves, D. B. Scariot, R. R. Guimaraes, C. V. Nakamura, R. S. Mendes and H. V. Ribeiro, PLoS One, 2016, 11, e0152092.
25. I. Munguira, I. Casuso, H. Takahashi, F. Rico, A. Miyagi, M. Chami and S. Scheuring, ACS Nano, 2016, 10, 2584.
26. S. Hanot, S. Lyonnard and S. Mossa, Nanoscale, 2016, 8, 3314.
27. D. Wang, C. He, M. P. Stoykovich and D. K. Schwartz, ACS Nano, 2015, 9, 1656.
28. M. J. Skaug, J. Mabry and D. K. Schwartz, Phys. Rev. Lett., 2013, 110, 256101.
29. G. Kwon, B. J. Sung and A. Yethiraj, J. Phys. Chem. B, 2014, 118, 8128.
30. E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield and D. A. Weitz, Science, 2000, 287, 627.
31. P. Charbonneau, Y. Jin, G. Parisi and F. Zamponi, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 15025.
32. M. J. Skaug, L. Wang, Y. Ding and D. K. Schwartz, ACS Nano, 2015, 9, 2148; M. J. Skaug and D. K. Schwartz, Ind. Eng. Chem. Res., 2015, 54, 4414.
33. B. Berkowitz, A. Cortis, M. Dentz and H. Scher, Rev. Geophys., 2006, 44, RG2003.
34. B. Wang, S. M. Anthony, S. C. Bae and S. Granick, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15160.
35. B. Wang, J. Kuo, S. C. Bae and S. Granick, Nat. Mater., 2012, 11, 481.
36. J. Guan, B. Wang and S. Granick, ACS Nano, 2014, 8, 3331.
37. K. He, F. B. Khorasani, S. T. Retterer, D. K. Thomas, J. C. Conrad and R. Krishnamoorti, ACS Nano, 2013, 7, 5122.
38. E. W. Montroll and G. H. Weiss, J. Math. Phys., 1965, 6, 167; H. Scher and E. W. Montroll, Phys. Rev. B: Solid State, 1975, 12, 2455.
39. J.-H. Jeon, E. Barkai and R. Metzler, J. Chem. Phys., 2013, 139, 121916.
40. T. Akimoto and E. Barkai, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 032915; T. Akimoto, Phys. Rev. Lett., 2012, 108, 164101; T. Geisel and S. Thomae, Phys. Rev. Lett., 1984, 52, 1936; T. Geisel, J. Nierwetberg and A. Zacherl, Phys. Rev. Lett., 1985, 54, 616; M. W. Deem and D. Chandler, J. Stat. Phys., 1994, 76, 911; C. Monthus and J.-P. Bouchaud, J. Phys. A: Math. Gen., 1996, 29, 3847; E. Bertin and J.-P. Bouchaud, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 026128; S. Burov and E. Barkai, Phys. Rev. Lett., 2007, 98, 250601; M. Dentz, et al. , Adv. Water Resour., 2012, 49, 13; T. Akimoto, E. Barkai and K. Saito, 2016, arXiv:1604.06175.
41. M. M. Meerschaert, E. Nane and Y. Xiao, Stat. Probab. Lett., 2009, 79, 1194; A. V. Chechkin, M. Hofmann and I. M. Sokolov, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 031112; J. H. P. Schulz, A. V. Chechkin and R. Metzler, J. Phys. A: Math. Theor., 2013, 46, 475001.
42. V. Tejedor and R. Metzler, J. Phys. A: Math. Theor., 2010, 43, 082002; M. Magdziarz, R. Metzler, W. Szczotka and P. Zebrowski, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 051103.
43. B. B. Mandelbrot and J. W. van Ness, SIAM Rev., 1968, 10, 422.
44. R. Kubo, Rep. Prog. Phys., 1966, 29, 255; P. Hänggi, Z. Phys. B: Condens. Matter Quanta, 1978, 31, 407; P. Hänggi and F. Mojtabai, Phys. Rev. E, 1982, 26, 1168; S. C. Kou, Ann. Appl. Stat., 2008, 2, 501; I. Goychuk, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 046125; I. Goychuk, Adv. Chem. Phys., 2012, 150, 187.
45. G. Pagnini, D. Molina-Garcia, T. M. Pham, C. Manzo and P. Paradisi, 2015, E-print arXiv:1508.01361.
46. W. Deng and E. Barkai, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 011112.
47. J. W. Haus and K. W. Kehr, Phys. Rep., 1987, 150, 263; S. Havlin and D. Ben-Avraham, Adv. Phys., 1987, 36, 695; G. C. Papanicolaou, “Diffusion in random media”, Surveys in applied mathematics, Plenum Press, New York, 1995, pp. 205–253; G. M. Zaslavsky, Phys. Rep., 2002, 371, 461; J. P. Bouchaud, A. Comtet, A. Georges and P. le Doussal, Ann. Phys., 1990, 201, 285.
48. A. J. Bray, S. N. Majumdar and G. Schehr, Adv. Phys., 2013, 62, 225.
49. Y. Mardoukhi, J.-H. Jeon and R. Metzler, Phys. Chem. Chem. Phys., 2015, 17, 30134.
50. P. Massignan, et al. , Phys. Rev. Lett., 2014, 112, 150603.
51. M. V. Chubynsky and G. W. Slater, Phys. Rev. Lett., 2014, 113, 098302.
52. T. Uneyama, T. Miyaguchi and T. Akimoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 032140.
53. R. Jain and K. L. Sebastian, J. Phys. Chem. B, 2016, 120, 3988.
54. J. Bewerunge, I. Ladadwa, F. Platten, C. Zunke, A. Heuer and S. Egelhaaf, Phys. Chem. Chem. Phys., 2016, 18, 18887.
55. W. Gotze and L. Sjogren, Rep. Prog. Phys., 1992, 55, 241.
56. R. Yamamoto and A. Onuki, Phys. Rev. Lett., 1998, 81, 4915.
57. R. Richert, J. Phys.: Condens. Matter, 2002, 14, R703.
58. L. Berthier and G. Biroli, Rev. Mod. Phys., 2011, 83, 587.
59. M. J. Saxton, Biophys. J., 1997, 72, 1744.
60. M. Platani, I. Goldberg, A. I. Lamond and J. R. Swedlow, Nat. Cell Biol., 2002, 4, 502.
61. Y. M. Wang, R. H. Austin and E. C. Cox, Phys. Rev. Lett., 2006, 97, 048302.
62. I. Goychuk and V. O. Kharchenko, Phys. Rev. Lett., 2014, 113, 100601.
63. M. Bauer, E. S. Rasmussen, M. A. Lomholt and R. Metzler, Sci. Rep., 2015, 5, 10072.
64. S. Hapca, J. W. Crawford and L. M. Young, J. R. Soc., Interface, 2009, 6, 111.
65. K. Yamamura, Popul. Ecol., 2002, 44, 93.
66. J. S. Clark, Am. Nat., 1998, 152, 204.
67. R. Nathan, Science, 2006, 313, 786.
68. K. Yamamura, Popul. Ecol., 2004, 46, 87.
69. G. T. Skalski and J. F. Gilliam, Ecology, 2000, 81, 1685.
70. J. K. M. Brown and M. S. Hovmoller, Sci. Transl. Med., 2002, 297, 537.
71. T. Kühn, T. O. Ihalainen, J. Hyvaluoma, N. Dross, S. F. Willman, J. Langowski, M. Vihinen-Ranta and J. Timonen, PLoS One, 2011, 6, e22962.
72. B. P. English, V. Hauryliuk, A. Sanamrad, S. Tankov, N. H. Dekker and J. Elf, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, E365.
73. C. B. Mast, S. Schink, U. Gerland and D. Braun, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 8030.
74. A. G. Cherstvy, A. V. Chechkin and R. Metzler, New J. Phys., 2013, 15, 083039; A. G. Cherstvy and R. Metzler, Phys. Chem. Chem. Phys., 2013, 15, 20220.
75. A. G. Cherstvy, A. V. Chechkin and R. Metzler, Soft Matter, 2014, 10, 1591.
76. A. G. Cherstvy, A. V. Chechkin and R. Metzler, J. Phys. A: Math. Theor., 2014, 47, 485002.
77. K. S. Fa and E. K. Lenzi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 061105.
78. K. S. Fa, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 020101(R).
79. M. Burgis, et al. , New J. Phys., 2011, 13, 043031.
80. M. Heidernätsch, On the diffusion in inhomogeneous systems, PhD thesis, TU Chemnitz, 2015.
81. H. Safdari, A. G. Cherstvy, A. V. Chechkin, F. Thiel, I. M. Sokolov and R. Metzler, J. Phys. A: Math. Theor., 2015, 48, 375002.
82. E. Geneston, R. Tuladhar, M. T. Beig, M. Bologna and P. Grigolini, 2015, arXiv:1601.02879.
83. S. C. Lim and S. V. Muniandy, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 021114.
84. R. Friedrich, J. Peinke, M. Sahimi and M. R. R. Tabar, Phys. Rep., 2011, 506, 87.
85. A. Bodrova, A. V. Chechkin, A. G. Cherstvy and R. Metzler, Phys. Chem. Chem. Phys., 2015, 17, 21791.
86. F. Thiel and I. M. Sokolov, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 012115.
87. J.-H. Jeon, A. V. Chechkin and R. Metzler, Phys. Chem. Chem. Phys., 2014, 16, 15811.
88. G. K. Batchelor, Math. Proc. Cambridge Philos. Soc., 1952, 48, 345.
89. A. G. Cherstvy and R. Metzler, J. Chem. Phys., 2015, 142, 144105.
90. A. Fulinski, J. Chem. Phys., 2013, 138, 021101; A. Fulinski, Acta Phys. Pol., 2013, 44, 1137.
91. A. Fulinski, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 061140.
92. A. G. Cherstvy and R. Metzler, J. Stat. Mech.: Theory Exp., 2015, P05010.
93. H. Safdari, A. V. Chechkin, G. R. Jafari and R. Metzler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 042107.
94. A. G. Cherstvy and R. Metzler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 012134.
95. A. W. C. Lau and T. C. Lubensky, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 011123.
96. A. Bodrova, A. V. Chechkin, A. G. Cherstvy and R. Metzler, New J. Phys., 2015, 17, 063038.
97. G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev., 1930, 36, 823.
98. A. Bodrova, A. V. Chechkin, A. G. Cherstvy, H. Safdari, I. M. Sokolov and R. Metzler, Sci. Rep., 2016, 6, 30520.
99. B. Lindner and I. M. Sokolov, Phys. Rev. E, 2016, 93, 042106.
100. C. Di Rienzo, V. Piazza, E. Gratton, F. Beltram and F. Cardarelli, Nat. Commun., 2014, 5, 5891.
101. S. Regev, N. Gronbech-Jensen and O. Farago, 2016, arXiv:1606.08632.
102. O. Farago and N. Gronbech-Jensen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 013301.
103. O. Farago and N. Gronbech-Jensen, J. Stat. Phys., 2014, 156, 1093.
104. J.-P. Bouchaud, J. Phys. I, 1992, 2, 1705.
105. G. Bel and E. Barkai, Phys. Rev. Lett., 2005, 94, 240602; A. Rebenshtok and E. Barkai, Phys. Rev. Lett., 2007, 99, 210601.
106. A. Lubelski, I. M. Sokolov and J. Klafter, Phys. Rev. Lett., 2008, 100, 250602.
107. S. Burov, R. Metzler and E. Barkai, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 13228.
108. Y. He, S. Burov, R. Metzler and E. Barkai, Phys. Rev. Lett., 2008, 101, 058101.
109. J.-H. Jeon and R. Metzler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 021103.
110. J.-H. Jeon, N. Leijnse, L. B. Oddershede and R. Metzler, New J. Phys., 2013, 15, 045011.
111. M. Schwarzl, A. Godec, E. Barkai and R. Metzler, unpublished.
112. A. Andreanov and D. S. Grebenkov, J. Stat. Mech.: Theory Exp., 2012, P07001.
113. J. Kursawe, J. Schulz and R. Metzler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062124.
114. J. H. P. Schulz, E. Barkai and R. Metzler, Phys. Rev. Lett., 2013, 110, 020602; J. H. P. Schulz, E. Barkai and R. Metzler, Phys. Rev. X, 2014, 4, 011028.
115. A. Godec and R. Metzler, Phys. Rev. Lett., 2013, 110, 020603.
116. J.-H. Jeon and R. Metzler, J. Phys. A: Math. Theor., 2010, 43, 252001.
117. P. Siegle, I. Goychuk and P. Hänggi, Phys. Rev. Lett., 2010, 105, 100602.
118. B. Lindner, New J. Phys., 2007, 9, 136; B. Lindner, J. Stat. Phys., 2008, 130, 523; B. Lindner, New J. Phys., 2010, 12, 063026.
119. P. Romanczuk, M. Bär, W. Ebeling, B. Lindner and L. Schimansky-Geier, Eur. Phys. J.: Spec. Top., 2012, 202, 1.
120. M. J. Saxton, Biophys. J., 2001, 81, 2226.
121. T. J. Feder, I. Brust-Mascher, J. P. Slattery, B. Baird and W. W. Webb, Biophys. J., 1996, 70, 2767.
122. P. Schwille, J. Korlach and W. W. Webb, Cytometry, 1999, 36, 176.
123. M. Weiss, H. Hashimoto and T. Nilsson, Biophys. J., 2003, 84, 4043.
124. A. Kusumi, Y. Sako and M. Yamamoto, Biophys. J., 1993, 65, 2021.
125. S. Hapca, J. W. Crawford, K. MacMillan, M. J. Wilson and I. M. Young, J. Theor. Biol., 2007, 248, 212.
126. C. Charalambous, G. Munoz-Gil, A. Celi, M. F. Garcia-Parajo, M. Lewenstein, C. Manzo and M. A. Garca-March, 2016, arXiv:1607.00189.
127. A. G. Cherstvy, A. V. Chechkin and R. Metzler, unpublished.

---