Communicated by Robert J. Silbey, Massachusetts Institute of Technology, Cambridge, MA, December 22, 2007 (received for review November 14, 2007)
Subdiffusive motion of tracer particles in complex crowded environments, such as biological cells, has been shown to be widespread. This deviation from Brownian motion is usually characterized by a sublinear time dependence of the mean square displacement (MSD). However, subdiffusive behavior can stem from different microscopic scenarios that cannot be identified solely by the MSD data. In this article we present a theoretical framework that permits the analytical calculation of first-passage observables (mean first-passage times, splitting probabilities, and occupation times distributions) in disordered media in any dimensions. This analysis is applied to two representative microscopic models of subdiffusion: continuous-time random walks with heavy tailed waiting times and diffusion on fractals. Our results show that first-passage observables provide tools to unambiguously discriminate between the two possible microscopic scenarios of subdiffusion. Moreover, we suggest experiments based on first-passage observables that could help in determining the origin of subdiffusion in complex media, such as living cells, and discuss the implications of anomalous transport to reaction kinetics in cells.
In the past few years, subdiffusion has been observed in an increasing number of systems (1, 2), ranging from physics (3, 4) or geophysics (5) to biology (6, 7). In particular, living cells provide striking examples for systems where subdiffusion has been repeatedly observed experimentally, either in the cytoplasm (6–9), the nucleus (10, 11), or the plasmic membrane (12–14). However, the microscopic origin of subdiffusion in cells is still debated, even if believed to be caused by crowding effects in a wide sense, as indicated by in vitro experiments (15–18).
The subdiffusive behavior significantly deviates from the usual Gaussian solution of the simple diffusion equation, and is usually characterized by a mean square displacement (MSD) (1) that scales as 〈Δr 2〉 ∼ t β with β < 1. Such a scaling law can be obtained from a few models based on different underlying microscopic mechanisms. Here, we focus on two possibilities (a third classical model of subdiffusion is given by the fractional Brownian motion that concerns processes with long-range correlations): (i) The first class of models that we consider stems from continuous time random walks (CTRWs) (1, 19) and their continuous limit described by fractional diffusion equations (1, 20). The anomalous behavior in these models originates from a heavy-tailed distribution of waiting times (21): at each step the walker lands on a trap, where it can be trapped for extended periods of time. When dealing with a tracer particle, traps can be out-of-equilibrium chemical binding configurations (22, 23), and the waiting times are then the dissociation times; traps can also be realized by the free cages around the tracer in a hard sphere-like crowded environment, and the waiting times are the life times of the cages (see Fig. 1 a). (ii) Another kind of model for subdiffusion relies on spatial inhomogeneities as exemplified by diffusion in deterministic or random fractals such as critical percolation clusters (24–26). The anomalous behavior is in this case caused by the presence of fixed obstacles (27) that create numerous dead ends, as illustrated by De Gennes's “ant in a labyrinth” (28) (see Fig. 1 b). These two scenarios can be classified as dynamic (CTRW) and static (fractal) in the nature of the underlying environment.
Although these two models lead to similar scaling laws for the MSDs, their microscopic origins are intrinsically different and lead to notable differences in other transport properties. This has strong implications, in particular, on transport-limited reactions (29), which will prove to have very different kinetics in the two situations. Because most functions of a living cell are regulated by coordinated chemical reactions that involve low concentrations of reactants [such as transcription factors or vesicles carrying targeted proteins (30)], and that are limited by transport, understanding the origin of anomalous transport in cells and its impact on reaction kinetics is an important issue.
Here, we describe and analytically calculate the following transport-related observables, based on first-passage properties, which allows, as shown below, discrimination between the CTRW and fractal models and permits a quantitative analysis of the kinetics of transport-limited reactions:
The first-passage time (FPT), which is the time needed for a particle starting from site S to reach a target T for the first time. This quantity is fundamental in the study of transport-limited reactions (31–33), because it gives the reaction time in the limit of perfect reaction. This quantity is also useful in target search problems (34–39) and other physical systems (40–42). We will be interested in both the probability density function (PDF) of the FPT, and its first moment, the mean FPT (MFPT).
The first-passage splitting probability, which is the probability to reach a target T 1 before reaching another target T 2, in the case where several targets are available. This quantity permits the study of competitive reactions (31).
The occupation time before reaction, which is the time spent by a particle at a given site T 1 before reaction with a target T 2. This quantity is useful in the context of reactions occurring with a finite probability per unit of time (43–45). We stress that the occupation time provides a finer information on the trajectory of the particle. In particular, the FPT is given by the sum over all sites of the occupation time. We will be interested in both the entire PDF of the occupation time and the mean occupation time.
On the theoretical level, our approach permits the direct evaluation of nontrivial first-passage characteristics of transport in disordered media in any dimension, whereas, so far, mainly effective one-dimensional geometries have been investigated (42). In particular, we calculate here MFPT, splitting probabilities and occupation time distribution of a random walk on percolation clusters, and discuss the potential implications of these results on reaction kinetics in living cells. These findings could lead to an experimental probing of the microscopic origin of subdiffusion in complex media such as living cells.
The article is organized as follows. In the first section, we set the theoretical framework and give explicit analytical expressions of the first-passage observables, which are summarized in Eqs. 13 – 15 . We then apply these results to the two above mentioned models of subdiffusion, namely the diffusion on fractal and CTRW models. In the second section, we discuss the relevance of these two models to describe anomalous transport in complex media such as living cells, and suggest experiments that could help in discriminating the microscopic origin of subdiffusion.
By using recent techniques developed in refs. 40, 46, and 47, we derive general analytical expressions of the first-passage observables. We consider a Markovian random walker moving in a bounded domain of N sites, with reflecting walls. Let W(r, t|r′) be the propagator, i.e., the probability density to be at site r at time t, starting from the site r′ at time 0, whose evolution is described by a master equation (48):
with a given transition operator ℒ. We denote by P(r,t |r′) the probability density that the first-passage time to reach r, starting from r′, is t. For the sake of simplicity we assume that the walker performs symmetric jumps and that the stationary distribution is homogeneous lim
t→∞
W(r, t|r′) = 1/N. The propagator and first-passage time densities are known to be related through (49)
Following ref. 40, this equation gives an exact expression for the MFPT, provided it is finite:
where
H
is the pseudo-Green function (50) of the domain:
It is also possible to compute splitting probabilities within this framework. If the random walker can be absorbed either by a target T
1 at r
1, or a target T
2 at r
2, a similar calculation yields:
where P
1 (resp. P
2) is the splitting probability to hit T
1 (resp. T
2) before T
2 (resp. T
1), and 〈T〉 is the mean time needed to hit any of the targets. This equation together with the similar equation obtained by inverting 1 and 2 and the condition P
1 + P
2 = 1, give a linear system of three equations for the three unknowns P
1, P
2, and 〈T〉, which can therefore be straightforwardly determined. In particular, the splitting probability P
1 reads:
where we used the notation H
ij = H(r
i|r
j). This formula extends a previous result (46, 47) obtained for simple random walks in the case of general Markov processes.
Beyond their own interest, the splitting probabilities allow us to obtain the entire distribution of the occupation time (45) N
i at site
i
for general Markov processes. Denoting P
ij(i|S) the splitting probability to reach i before j, starting from S, we have P(N
i = 0) = P
iT(T|S), and for k ≥ 1:
where
and E
2 is the probability to reach T starting from i without ever returning to i, which reads (45):
In particular, the mean occupation time is then given by
We stress that Eq. 7
gives the exact distribution of the occupation time for all regimes. It follows in particular that the large time asymptotics of the occupation time distribution is exponential. Actually one can argue in the general case that the FPT is also exponentially distributed at long times. This comes from the fact that the transition operator ℒ has a strictly negative discrete spectrum for a finite volume N (see ref. 48).
Eqs. 3
,
6
, and
10
give exact expressions of the first-passage observables as functions of the pseudo-Green function H. The key point is that, as shown in ref. 40, H can be satisfactorily approximated by its infinite space limit, which is precisely the usual Green function G
0:
where W
0 is the infinite space propagator. Following ref. 40, we assume that the problem is scale invariant and we use for W
0 the standard scaling (24):
where the fractal dimension d
f characterizes the accessible volume V
r ∼ rdf
within a sphere of radius r, and the walk dimension dw
characterizes the distance r ∼ t1/dw
covered by a random walker in a given time t. The form (Eq. 12
) ensures the normalization of W
0 by integration over the whole fractal set. Note that the MSD is then given by 〈Δr
2〉 ∼ t
β with β = 2/dw
. A derivation given in ref. 40 then allows us to extract the scaling of the pseudo-Green function H, and eventually yields for the MFPT:
where explicit expressions of A and B are given in ref. 40. We stress that, in the case of compact exploration (dw
> df
), the MFPT depends on a single constant B. Indeed, the constant A introduced in ref. 40 can be shown to be actually 0 in this case of compact exploration in scale invariant media. In fact, the previous analysis of the pseudo-Green functions also permits us to obtain explicit expressions of the splitting probabilities and mean occupation times:
and
where r
ij
= |r
i
− r
j
| is different from 0. Note that the entire distribution of Ni
is obtained similarly by estimating E
1 and E
2 as defined by Eqs. 8
and
9
. Strikingly, the constants A and B do not depend on the confining domain and can be written solely in terms of the infinite space scaling function ∏. We point out that in the case of compact exploration the expression of the splitting probability is fully explicit and does not depend on ∏. Eqs. 13
–
15
therefore elucidate the dependence of the first-passage observables on the geometric parameters of the problem, and constitute the central theoretical result of this article. We discuss the implications of these results on explicit examples in the next paragraph.
Critical percolation clusters (see Fig. 1 b) constitute a representative example of random fractals (24, 25, 51). Here, we consider the case of bond percolation, where the bonds connecting the sites of a regular lattice of the d–dimensional space are present with probability p. The ensemble of points connected by bonds is called a cluster. If p is above the percolation threshold p c, an infinite cluster exists. If p = p c, this infinite cluster is a random fractal characterized by its fractal dimension df . We consider a nearest-neighbor random walk on such a critical percolation cluster, with the so-called “blind ant (49)” dynamics: on arrival at a given site s, the walker attempts to move to one of the adjacent sites on the original lattice with equal probability. If the link corresponding to this move does not exist, the walker remains at site s. This walk is characterized by the walk dimension dw . In the example of the three-dimensional cubic lattice, one has d f = 2.58 … , and dw = 3.88… (25) and the motion is subdiffusive with β = 2/d w ≃ 0.51. For a given critical percolation cluster, namely, for a given configuration of the disorder, the theoretical development of previous paragraph holds, and the first-passage observables are given by the exact expressions (Eqs. 3 , 6 , and 10 ). However, the variations between different realizations of the disorder have to be taken into account, and averaging has to be performed to obtain meaningful quantities. It is shown in Materials and Methods that expressions 3, 6, and 10 actually still hold after disorder averaging.
Two scenarios of subdiffusion for a tracer particle in crowded environments. (a) Random walk in a dynamic crowded environment. The tracer particle evolves in a cage whose typical life time diverges with density. This situation can be modeled by a CTRW with power-law distributed waiting times. (b) Random walk with static obstacles. This situation can be modeled by a random walk on a percolation cluster.
Fig. 2 shows that the simulations fit very well the expected scaling. Both the volume dependence and the source–target distance dependence are faithfully reproduced by our theoretical expressions, as shown by the data collapse of the numerical simulations.
Numerical simulation of first-passage observables for random walks on three-dimensional percolation clusters. All of the embedding domains have reflecting boundary conditions. (a) MFPT for random walks on 3-dimensional critical percolation clusters. For each size of the confining domain, the MFPT, normalized by the number of sites N, is averaged both over the different target and starting points separated by the corresponding chemical distance, and over percolation clusters. The black plain curve corresponds to the prediction of Eq. 15 with dw c − dc f ≃ 1. (b) Splitting probability for random walks on 3-dimensional critical percolation clusters. The splitting probability P 1 to reach the target T 1 before the target T 2 is averaged both over the different target points T 2 and over the percolation clusters. The chemical distance ST 1 = 10 is fixed whereas the chemical distance ST 2 = T 1 T 2 is varied. The black plain curve corresponds to the explicit theoretical expression 14 with d w c − df c ≃ 1. (c) Occupation time for random walks on critical percolation clusters. For each size of confining domain, the occupation time of site T 1 before the target T 2 is reached for the first time is averaged over the different target points T 2 and over the percolation clusters. The chemical distance ST 1 = 10 is fixed whereas the chemical distance ST 2 = T 1 T 2 is varied. The black plain curve corresponds to the prediction of Eq. 15 with dw c − df c ≃ 1. (d) The MFPT for random walks on percolation clusters above criticality for a 25 × 25 × 25 confining domain. The MFPT, normalized by the number of sites N, is averaged both over the different target and starting points separated by the corresponding chemical distance, and over the percolation clusters.
If the bond concentration p is above the percolation threshold p c, a correlation length ξ∝(p − p c) −ν appears, where ν = 0.87 for d = 3. At length scales smaller than ξ, the percolation cluster is fractal, with the same fractal dimension df as the critical percolation cluster, and diffusion is anomalous. At length scales larger than ξ, the fractal dimension of the percolation cluster recovers the space dimension d and diffusion is normal (24).
Along the lines of the previous section, we thus expect the pseudo-Green function H to scale as r
dw
−df
for r < ξ, and as r
d−2 for r > ξ. More explicitly, on the example of the MFPT, we expect for the three-dimensional cubic lattice
Similarly, the other first-passage observables display a cross-over between these two regimes around ξ. The simulations do show very well the transition between the two regimes (see Fig. 2
d).
The CTRW is not necessarily Markovian, unlike the fractal case, and therefore the methodology above cannot be applied directly. The distribution of the FPT for CTRWs was however obtained recently in ref. 52. We here briefly recall these results, and derive analytical expressions of the other observables. The CTRW is a standard random walk with random waiting times, drawn from a PDF ψ(t). The CTRW model has a normal diffusive behavior if the mean waiting time is finite. For heavy-tailed distributions such that
the mean waiting time diverges for α < 1 and the walk is subdiffusive since the MSD scales like 〈Δr
2〉 ∼ t
β with β = α (see refs. 1 and 3). Here, τ is a characteristic time in the process. We focus on the representative case of a one-sided Levy stable distribution (49) ψ(t), which satisfies Eq. 17
and whose Laplace transform is ψ̂(u) = exp(−τ
α
u
α
) (0 < α < 1).
We now derive the relation between the FPT to the site r
T
, starting from r
S
for the standard discrete-time random walk and the CTRW. Denoting π(t) the probability density of the FPT for the CTRW, and Q(n) the probability density of the FPT for the discrete-time random walk, n being the number of steps, one has
which is conveniently rewritten after Laplace transformation as
where Q̂(z) = Σn=1
∞
Q(n)zn
is the generating function of the discrete-time random walk.
Several comments are in order. (i) The small u limit shows that the MFPT is infinite, and the long-time behavior of π(t) is directly related to the MFPT of the discrete-time simple random walk:
It should be noted that as soon as Q̂(z) is exactly known (such as for d = 3 in the large N limit; see ref. 52), the entire distribution of the FPT can be obtained. (ii) Because splitting probabilities are time-independent quantities, they are exactly identical for CTRW and standard discrete-time random walks, and are therefore given by Eq. 14
with the space dimension d and the walk dimension dw
= 2. (iii) The same decomposition as Eqs. 18
and
19
holds for the distribution π
i(t
i) of the occupation time t
i of site i, where the distribution of the occupation time F(N
i) for the discrete-time random walk has to be introduced. This yields
Interestingly, as F(N
i
) is explicitly given by Eq. 7
, the entire distribution of the occupation time can be derived.
We emphasize that a proper definition of the mean values of the first-passage observables (namely the MFPT and the mean occupation time) is provided by introducing a truncated distribution (with cutoff t c) of waiting times in place of ψ(t). Because this allows the definition of a mean waiting time τm = C ∫0 tc tψ(t)dt (where C normalizes the truncated PDF), the MFPT is then given by 〈T〉 = τ m 〈n〉, and the mean occupation time reads 〈t i〉 = τ m 〈N i 〉.
Note that our results show that the first-passage observables scale with the geometric parameters N and r exactly as a simple random walk. Their scaling dependence is therefore given by Eqs. 13 – 15 , where df is the space dimension d and dw = 2.
We first discuss the relevance of the two models, CTRW and diffusion, on fractals to describe anomalous transport in confined systems such as the cytoplasm and membrane of living cells. The cell is known to be a highly complex and inhomogeneous molecular assembly, composed of numerous constituents that may vary widely from one cell type to another. Here, we wish to distinguish between two types of effects on transport in cellular medium. First, the overall density of free proteins and molecular aggregates is very high, be it in the cytoplasm or in the plasma membrane. In such a crowded environment, a tracer particle is trapped in dynamic “cages” whose life times are broadly distributed at high densities and leading to Eq. 17 . This dynamic picture therefore fits the hypothesis of the CTRW model. Second, the cytoskeleton is made of semiflexible polymeric filaments (such as F-actin or microtubules) that can be branched and cross-linked by proteins. This scaffold therefore acts as fixed obstacle constraining the motion of the tracer. Moreover, the cytoplasm can be compartmentalized by lipid membranes that further constrain the tracer. Such environment with obstacles can be described in a first approximation by a static percolation cluster. How could one discriminate between these two mechanisms having markedly different physical origins?
The first-passage observables derived earlier make it possible to distinguish between the two models of subdiffusion, as summarized in Table 1. (i) The first-passage time has a finite mean and exponential tail for the fractal model, whereas it has an infinite mean and a power-law tail in a CTRW model. Analyzing the tail of the distribution of the FPT therefore provides a first tool to distinguish the two models. Because experiments can only find the first-passage up to a certain time, we need to use the above-mentioned truncated means to define the MFPT for CTRW. In this case, the scaling of the MFPT for CTRW with the source–target distance is the same as for a simple random walk, and can be distinguished from the scaling of the MFPT on random fractals. These two scalings are strikingly different for d = 3: the CTRW performs a noncompact exploration of space (dw = 2 < 3 = d) leading to a finite limit of the MFPT at large source–target distance, whereas exploration is compact for a random walker on the percolation cluster (dw > df ) leading to a scaling ∝ rdw−df of the MFPT. We highlight that this feature could have very strong implications on reaction kinetics in cells. Indeed, in the cases where the fractal description of the cell environment is relevant, our results show that reaction times crucially depend on the source–target distance r. The biological importance of such dependence on the starting point was recently emphasized in ref. 39 on the example of gene colocalization. However, when the CTRW description of transport is valid, reaction times do not depend on the starting point at large distance r. (ii) The splitting probabilities for the CTRW model and for the fractal models have different scalings with the distance between the source and the targets. As mentioned previously the difference is more pronounced for d = 3: the probability to reach the furthest target T 2 vanishes as r−(dw−df) for the fractal model, r being the distance ST 1 with the notations of Fig. 3, but it tends to a constant for d = 3 according to the CTRW model. As discussed earlier, this could have important consequences for the kinetics of competitive reactions in cells. (iii) As for the occupation time, both its distribution and the scaling of the conditional mean with the distances ST 1 and ST 2 can be used to distinguish between models. The advantage of the mean occupation time is that it can still discriminate between the models after averaging over initial conditions, and could therefore be used even with a concentration of tracers.
Comparison of first-passage observables for CTRW and fractal models for d = 3
Schematic proposed set-up to measure first-passage observables.
We now briefly discuss potential experimental utilizations of first-passage observables. The schematic setup that we propose to measure these observables relies on single-particle tracking techniques (see Fig. 3). We consider a single tracer, either a fluorescent particle or a nanocrystal, moving in a finite volume such as a living cell, a microfluidic chamber, or vesicle. A laser excitation defines the starting zone S. As soon as the tracer enters S a signal is detected and a clock is started. Similarly, a second laser excitation defines the target zone T 1 and allows the measurement of the FPT of the tracer at T 1. In the same way, a third laser can detect a second target T 2: counting the time spent by the tracer in T 2 before reaching T 1 gives exactly the occupation time. Splitting probabilities are straightforwardly deduced.
Finally, this theoretical framework can be extended to cover more realistic situations. First, subdiffusion could result in some systems from a combination of both the dynamic (CTRW) and static (diffusion on fractal) mechanisms. Interestingly, our approach can be adapted to study the example of CTRWs on a fractal that models such situations (54). Indeed, the same decomposition as in Eq. 18 holds in this case and shows that the dependence of the first-passage observables (defined with truncated means if needed) on the source–target distance is exactly the same as in the case of a standard discrete-time random walk on the fractal, and therefore gives access to the dimensions dw and df of the fractal. In turn, the tail of the distribution of the FPT is in this case reminiscent of the single-step waiting time distribution defining the CTRW as shown by Eq. 20 (see also ref. 54). First-passage observables therefore permit us, in principle, to isolate and characterize each of the CTRW and fractal mechanisms even when they are both involved simultaneously. Second, in various systems subdiffusion occurs over a given time scale or length scale, crossing over to the regular diffusive behavior. Both models can be adapted to capture this effect. In the fractal model the fractal structure persists up to the cross-over length scale (which is the correlation length ξ in percolation clusters above criticality), and the waiting time distribution for the CTRW model has a Levy-like decay until the cross-over time scale, after which the decay is faster so that the mean waiting time becomes finite. The MFPT will exist in both of these modified models, but the CTRW model leads to a normal scaling of the MFPT with the volume and the source–target distance: namely, it corresponds to the results of the simple random walk, with the same time step as the mean waiting time. However, a truncated fractal structure would lead to the same scaling on larger scales, but to a scaling as in Eq. 15 at smaller scales. The small-distance behavior of the MFPT can thus discriminate the two models. The same conclusion holds for the splitting probabilities and occupation times: the small-length behavior will also differ.
Our approach therefore permits us to explore the scaling of first-passage observables for two representative models of subdiffusion as a methodology to discriminate between underlying mechanisms for subdiffusion and to gain insight into the microscopic origin of subdiffusion and the nature of transport-limited reactions in complex systems.
We will denote by X̄ the average of X over the disorder, and assume that all configurations have the same volume N, which is a nonrestrictive condition in the large N limit since N is self-averaging. Eqs. 3 , 6 , and 10 then show that averaging the first-passage observables amounts to averaging the pseudo-Green function, and therefore the propagator in virtue of Eq. 4 . In the case of a random walk on a critical percolation cluster it has been shown that the propagator has a multifractal behavior (25). This means that the propagator W(r, t) has a very broad distribution, and is not self-averaging: its typical value is not its average value, which is dominated by rare events. In particular, a scaling form of the averaged propagator is not available. However, this difficulty can be bypassed if one considers the chemical distance x, i.e., the step length of the shortest path between two points. Indeed, in the chemical space, the propagator does have a simple fractal scaling (25, 51) and in the infinite volume limit the averaged propagator W̄0(x, t) satisfies the scaling form 12 (see ref. 25). Note that this property is shared by most random fractals (25), and makes the chemical distance space a powerful tool to calculate disorder averages. The formalism derived in the previous section can therefore be used, and the scaling laws of the MFPT, splitting probability, and mean occupation time averaged over the disorder are given in chemical space by Eqs. 13 – 15 , where r is to be replaced by the chemical distance x. Note that, in the chemical space, the fractal dimension is given by df c = d f /d min and walk dimension is dw c = d w /d min . The dimension d min is the fractal dimension of chemical paths and permits us to recover the dependence on the Euclidian distance r through the scaling (24) x ∼ r dmin, with d min = 1.24 in the case of the three-dimensional cubic lattice (24).
These scaling laws for the first-passage observables can be tested numerically. We simulated in Fig. 2 several critical percolation clusters on the three-dimensional cubic lattice embedded in the confining domain, and we averaged for each set of chemical distances {xij } the desired observable over all configurations of source and targets yielding the same set {xij }.
We thank P. Desbiolles for useful discussions.
Author contributions: S.C., V.T., R.V., O.B., and J.K. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.
The authors declare no conflict of interest.
