Ergodic and non-ergodic anomalous diffusion in coupled stochastic processes

Golan Bel¹ and Ilya Nemenman

Center for Nonlinear Studies and Computer, Computation and Statistical Sciences Division, Los Alamos National Laboratory, Los Alamos NM 87545, USA

¹ Author to whom any correspondence should be addressed.

E-mail: golanbel@lanl.gov

Received 9 May 2009
Published 12 August 2009

Contents

• 1. Introduction
• 2. The model
• 3. Time averaged MSD (TAMSD)
• 4. Discussion
• Acknowledgments
• References

1. Introduction

Single-molecule kinetics has come within reach of biophysical experiments [1]–[3], and theoretical and computational tools for analysis of such processes have experienced a corresponding growth [4]–[8]. It is clear that the combinatorially large number of microscopic steps involved in even the simplest of biochemical events makes their rigorous stochastic treatment difficult. For example, gene expression, often modeled as a single-step mRNA creation, in fact, includes transcription-factor-DNA binding, polymerase recruitment, transcriptional bubble formation and multiple elongation steps, each of which is a complex process on its own.

Therefore, any theoretical analysis of stochastic biochemical processes necessarily involves coarse-graining: identifying a small subset of dynamical variables that are modeled explicitly, while agglomerating the rest into an effective behavior. Such coarse-grained dynamics is often modeled using the master equation or the Langevin approaches, which require Markoviness or white-noise random forcing. Both of these assumptions are, generally, flawed, and quantitative corrections have been worked out in certain cases [9, 10]. Less explored is the possibility that internal degrees of freedom can introduce qualitative differences, such as long-range temporal correlations among state transitions and non-diffusive behavior of the observed quantities.

A well-studied example shows that this is possible for a random walk on a discrete lattice. In [11], Weiss and Havlin analyzed a two-dimensional diffusion model, known as the comb model. There dynamics along the y-coordinate is unlimited, while motion along x is allowed only when y = 0. This results in x = 0 and , that is, in a subdiffusive motion of x. This model is hardly realistic in a biochemical context due to the discontinuous dependence of the diffusion coefficient on y. However, it is plausible that diffusive dynamics of a real biological or chemical variable in the state space or in the physical space depends on unobserved, decimated variables in some other nontrivial way. For example, in a chemotaxing E. coli, the number of unobserved signaling proteins CheY-P is coupled to the distribution of times a bacterial motor rotates counterclockwise, and the bacterium swims straight. For a fixed concentration of CheY-P, obtained by modifying the chemotaxis network, the distribution is essentially exponential [12], resulting in a regular diffusive motion of the bacterium. But in a wild-type bacterium, as the number of CheY-P fluctuates (diffuses in the number space even for a constant external signal), the distribution becomes a power law, and bacteria exhibit superdiffusive real-space motion. While not true in this particular system, the distribution of clockwise rotation times could have been strongly coupled to CheY-P as well. This would have resulted in a power-law distribution of times that the bacterium spends reorienting itself without moving forward, and hence in its subdiffusive motion. In both cases, neglecting the CheY-P fluctuations and describing bacterial motion as a normal diffusion is qualitatively wrong.

In this paper, we abstract out the detailed biology and explore these types of phenomena from the point of view of statistical physics. We derive the properties of a diffusion process, for which the diffusion coefficient depends on the state of another, unobserved, variable. We show that, quite generally, such dependence leads to anomalous diffusion of the observed process, suggesting that traditional stochastic approaches may fail, and that more thought should be given to modeling stochastic phenomena in complex interacting systems, in particular in biophysics.

2. The model

Our model is described by two variables x and y, which may represent, in particular, concentrations of two interacting chemical species. The x is considered as the observed quantity, while y is assumed hidden (i.e. unobserved). Particles of both species can be created and destroyed, which results in an overdamped diffusive motion of the system in the concentration space (we disregard the directional drift for simplicity, but it can be reintroduced easily). We assume that the diffusion of x is y-dependent. That is,

Here γ_x and γ_y are the effective friction coefficients (assumed to be homogeneous) corresponding to x, y; η(t) and ξ(t) are independent, zero-mean white noise forces such that

The idea of multiplicative noise was explored before. In [13], the authors considered the case of a Langevin process (not overdamped), in which a function of the velocity serves as a `filter' for the white noise. In [14], the multiplicative noise enters as a random friction coefficient; however, the distribution of the random friction is independent of time. Similarly in [15], the noise enters as a random mass with distribution which is again independent of time. Our model is different from the above mentioned models since the distribution of the random coupling parameter C(y) is time dependent; moreover it depends in an arbitrary way (through the function C(y)) on another variable, which diffuses and hence has long-range temporal correlations. Further, the coupling considered in our model does not introduce a directional bias since it is multiplied by a white noise ξ, which takes both positive and negative values.

The PDF of y is that of a normal diffusion

where the initial condition is y(t = 0) = y₀.

The dynamics of the mean square displacement (MSD) x(t)², where ··· stands for the average over the white noises η and ξ, can be derived. Formally integrating equation (2) and substituting it in the expression for the derivative of x² yields,

Averaging over the noise ξ(t) yields the dynamics of the MSD of x conditional on y(t),

To get the marginal expectation x(t)², we now average the conditional expectation over y, which is distributed as in equation (5):

The function C(y) may take different forms for different systems. The simplest case is when the dynamics of x is independent of y and C(y) = C = const. Substituting this into equation (8) yields the expected trivial result .

Another scenario is that x can evolve in time only for a given range of y values (|y| < y₁), which resembles the discrete comb model of Weiss and Havlin [11]. Notice, however, that the comb model is based on geometric constraints, i.e. the teeth, while in our case the coupling between the two process is not due to the topology, but rather due to the physical nature of the processes. The similarity arises due the fact that in both models the first passage time distribution in the infinitely long y-axis, dominate the dynamics. Indeed, substituting C(y(t)) = CΘ(y₁–y(t))Θ(y₁ + y(t)), where C is a dimensionless constant, and Θ(y) is the Heaviside Theta function, into equation (8), we see that, for t  y₁²/(12D_y), in agreement with [11]. If C(y) falls off exponentially as y→∞, the same subdiffusion exponent is recovered.

A more interesting case is when, at large y, C(y) falls off, but not too sharply. We consider a power-law form, namely

where A is a constant with the units of inverse length, and α is a dimensionless parameter. This form corresponds to repressive, cooperative Hill kinetics, which describes many biochemical processes [16]. Typical diffusive trajectories with this C(y), A = 1, and α = 0, 0.625 are shown in figure 1.

Assuming that the behavior of C(y) for large y (i.e. C(y)~|y|^–α) dominates the t→∞ dynamics of , a simple scaling argument suggests that . However, this is clearly wrong for large α, suggesting that must pick up an anomalous scaling due to the y→0 properties of C(y). In what follows, we derive the long-time behavior of in a more rigorous way.

Equation (8) with C(y) as in equation (9), and considering the long-time limit ( ) gives

In the first integral, we approximate the integrand as a constant for t→∞, and, in the second integral, we neglect 1 compared with , thus obtaining

where is the incomplete Gamma function. Integrating equation (11) over t results in the long-time behavior of

where D_{1, 2} are constants depending on the model parameters D_x, D_y, α and A. This implies that, for α < 1/2, the long-time behavior is dominated by , as the scaling argument suggests. However, for α > 1/2, the scaling argument breaks and . Note that when the C(y) falls faster than , the diffusion exponent is exactly the same as in the case in which the dynamics of x is limited to a finite range near y = 0.

The case of α = 1/2 is special and the integral of equation (10) can be calculated exactly, yielding

where G denotes the Meijer G function [17]. The leading order term of the MSD is then

The coupling function C(y) depends only on the variable y, which undergoes normal diffusion. Thus we can derive the probability distribution of C(y). In order to simplify the notation, in what follows we use C to denote the coupling function without explicitly specifying its dependence on y. When C is given by equation (9), its probability density is

where 0≤C≤1. In the above derivation, we set y(t = 0) = 0 for simplicity. We see that the distribution of the noise introduced by C is time dependent. Note that the temporal correlations of C are also important for the dynamic properties of x.

So far we have considered only situations in which the motion of x was slowed at large y, but we can also consider the opposite scenarios, when large y promotes diffusion in x, as in [12]

This now resembles the Hill activation kinetics for low concentration of the substrate molecules [16]. Typical trajectories for this form of the coupling function with β = 0.75 are shown in figure 1.

For coupling of this form, the dynamics of the MSD (equation (8)) yields

which, in the long-time limit, gives , and

The distribution of the coupling function C in this case is given by

where C≥0.

To confirm our analytical results, we performed numerical simulations for the different cases considered above. In figure 2, we present a comparison of the diffusion exponent ν (defined by ) versus the coupling parameter κ (for the subdiffusion scenario, equation (9), κ = α, and for the superdiffusion scenario, equation (16), κ = –β) between the simulations and the analytical results. The simulations were done according to equations (1) and (2) with γ_{x, y} = 1, D_{x, y} = 1 and dt = 1. We averaged the results over 10⁴ trajectories, each of duration 10⁷,  ..., 10⁸ dt.

Figure 2. Leading order diffusion exponent ν defined by for the coupled stochastic processes model with different couplings between the diffusion of x and y, measured by the exponent κ. For the case of hindered diffusion of x, equation (9), κ = α, while for the enhanced diffusion, equation (16), κ = –β. Numerical simulations (points) and theoretical predictions (line) agree for both scenarios.

A simple linear regression to was performed to estimate ν. Since the standard parameter errors obtained for the regressions were negligible, the error bars of ν were estimated from the variability of the fitted values as we changed the domain of t, for which the fits were performed. Figure 2 shows a clear agreement between our theoretical results and the simulations. Note that, in certain cases, the convergence to the leading behavior of as t→∞ is slow since the difference between the exponents of the leading and the subleading terms is small. This slowness determined the lengths of the simulations.

3. Time averaged MSD (TAMSD)

There are many models of anomalous diffusion, including a time dependent friction coefficient in the Langevin equation [18], continuous time random walk (CTRW) [19], fractional Brownian dynamics [20], fractional Langevin dynamics [21]–[23] and Langevin dynamics with colored noise [24], to name a few. For a new model resulting in an anomalous diffusion, it is important to see if it can be reduced to one of these more familiar constructions. For example, the t→∞ behavior of the original comb model is equivalent to CTRW [25] with a power-law tail of the distribution of the times between successive jumps along x. However, in our model, the analogy is not as straightforward since the continuous dynamics of y induces temporal correlations among successive motions along x.

To understand the relation of the coupled diffusion model to the others in the literature, we note that all of them yield the same behavior for the ensemble averaged MSD in the long-time limit. However, they still differ from each other in the short time behavior, the shape of the distribution, and even in the long time behavior of time-averaged quantities (for example, the CTRW exhibits ergodicity breaking [26]). In particular, the TAMSD is an important property (it is the TAMSD that is observed in typical single-molecule diffusion experiments in biological systems [3, 27], and the number of recorded trajectories is often insufficient to estimate ensemble averages).

The TAMSD is defined as

which averages the squared displacement of a particle in time Δ (the time lag) over a single-particle trajectory of duration t. For CTRW, the TAMSD is a random quantity and even its ensemble average still exhibits aging, that is, dependence on the measurement duration [28, 29] t in equation (20). On the contrary, for the fractional Brownian and Langevin dynamics, the TAMSD converges to the ensemble average MSD for long times [30].

We investigated the behavior of the TAMSD in our model with repressive coupling numerically. We find that, when the scaling argument holds, namely for α≤1/2 (see equation (9)), the TAMSD is not a random quantity, but it still shows aging, as we would expect for Langevin dynamics with a time-dependent friction. On the other hand, when α≥1/2, and the diffusion exponent is ν = 1/2, the TAMSD shows a similar behavior to that of the CTRW [28]. In figure 3, we show the TAMSD versus the time lag Δ for α = 0.75 and 0.25. Each line shows the TAMSD for a single trajectory. The solid red lines represent trajectories of duration t = 10⁵, and the dashed blue lines are for t = 10⁶. All lines show a linear dependence on Δ, but with different coefficients. For α = 0.25, the TAMSD lines converge as the trajectory duration grows (the dashed blue lines essentially overlap), while for α = 0.75, the lines remain random. This is a clear indication of ergodicity breaking in our model for α > 0.5. Further, this analysis suggests that the coupled Langevin processes model stands as its own class among other anomalous diffusion models, exhibiting time-dependent diffusion coefficient Langevin dynamics for certain forms of coupling, and some aspects of ergodicity-breaking CTRW for others.

Figure 3. The TAMSD versus Δ for α = 0.25, 0.75 (left and right panels, respectively). All other parameters are the same as for the previous figures. We used 20 trajectories of duration t = 105 (solid red lines) and 20 trajectories of duration t = 106 (dashed blue lines). The straight lines with slope 1 in the log-log scale reflect a linear dependence of the TAMSD on the time lag Δ. For α = 0.25, the TAMSD of different trajectories converges as the measurement time increases (the blue lines collapse, and the red do not), namely a time average of a single particle behaves like the ensemble average. However, for α = 0.75, there is no such convergence, just like for CTRW (the coefficient of proportionality between and Δ varies from trajectory to trajectory). In both cases, the ensemble average of the TAMSD decreases as the measurement time increases, indicating aging.

4. Discussion

In this paper, we introduced a model of coupled diffusive processes, where the diffusion of the observed variable x is coupled to the value of a hidden variable y. We showed that the dynamics of x exhibits anomalous diffusion for every considered form of the coupling between the variables. Depending on the nature of the coupling, the motion of x can be sub or superdiffusive (and even superballistic, as is the case of a frictionless particle subject to a white noise). Further, even for an arbitrary `strong' (such that the dynamics of x is limited to a small range of y values) repressive xy coupling, the x diffusion exponent ν is limited from below by 1/2 (anomalous scaling), so that localization of x is impossible. Even though the long-time ensemble-averaged behavior of our model is similar to that of many others describing anomalous diffusion, the model does not reduce to any of them, exhibiting an effective time-dependent diffusion coefficient, aging and ergodicity breaking for different values of its parameters.

The anomalous scaling and the ergodicity breaking appear for `strong' xy coupling (i.e. large α). This is because, for α < 1/2, motion of particles away from y = 0 contributes substantially to the ensemble-averaged MSD of x. On the contrary, for α > 1/2, only motion near y = 0 is important, namely the first passage time in an infinite one-dimensional system (the diffusion process in y) plays a major role in the dynamics of the observed quantity, x. It was also verified numerically that, for α > 1/2, where the dynamics of x is dominated by motion in a narrow range near y = 0, the propagator is non-Gaussian as expected from a CTRW-like renewal process. A similar result holds for the superdiffusive regime, β > 0. On the other hand, when the motion at any y is important (0≤α < 1/2), the process in x has long-time correlations and the propagator takes a Gaussian form (similar to that of a diffusion process with a time-dependent diffusion coefficient).

While important in its own right, the coupled diffusion model raises its most important questions in the biological domain. Unobserved dynamical quantities lead to anomalous diffusion in E. coli chemotaxis [12], or in mRNA diffusion in cells [3]. Further, some of the best established models of cellular regulation involve coarse-graining of dynamics. For example, in some models of the E. coli lac operon, the lac repressor itself is an unobserved variable [31], which is coupled to the speed of production of the lactose permease and the lactose-utilizing enzyme and, through them, to the import and the degradation of lactose in the cell. Since any coupling may lead to anomalous diffusion and in some cases even to ergodicity breaking, it begs the question whether relying on the common Langevin or master equation analysis of stochasticity of the lac repressor or other regulatory circuits, such as the λ-phage [4, 8], mar [1], and others, is rigorously justifiable.

Acknowledgments

We thank E Barkai and A Zilman for stimulating discussions and LANL Center for Nonlinear Studies for support and for striving to maintain an island of sanity in the increasingly more complex environment. This work was funded by LANL Laboratory Directed Research and Development program.

References

[1]

Guet C C et al 2008 Minimally invasive determination of mrna concentration in single living bacteria Nucleic Acids Res. 36 e73 
CrossRefPubMed

[2]

English B P et al 2006 Ever-fluctuating single enzyme molecules: Michaelis–Menten equation revisited Nat. Chem. Biol. 2 87 
CrossRefPubMed

[3]

Golding I and Cox E C 2006 Physical nature of bacterial cytoplasm Phys. Rev. Lett. 96 098102 
CrossRefPubMed

[4]

McAdams H H and Arkin A 1998 Simulation of prokaryotic genetic circuits Annu. Rev. Biophys. Biomol. Struct. 27 199–224 
CrossRefPubMed

[5]

Paulsson J 2004 Summing up the noise in gene networks Nature 427 415–8 
CrossRefPubMed

[6]

Sinitsyn N A, Hengartner N and Nemenman I 2009 Adiabatic coarse-graining and simulations of stochastic biochemical networks Proc. Natl Acad. Sci. USA 106 10546 
CrossRefPubMed

[7]

Sasai M and Wolynes P G 2003 Stochastic gene expression as a many-body problem Proc. Natl Acad. Sci. USA 100 2374–9 
CrossRefPubMed

[8]

Aurell E and Sneppen K 2002 Epigenetics as a first exit problem Phys. Rev. Lett. 88 048101 
CrossRefPubMed

[9]

Gopich I V and Szabo A 2006 Theory of the statistics of kinetic transitions with application to single-molecule enzyme catalysis J. Chem. Phys. 124 154712 
CrossRefPubMed

[10]

Sinitsyn N A and Nemenman I 2007 The berry phase and the pump flux in stochastic chemical kinetics Europhys. Lett. 77 58001 
IOPscience

[11]

Weiss G H and Havlin S 1986 Some properties of a random walk on a comb structure Physica A 134 474–82 
CrossRef

[12]

Korobkova E et al 2004 From molecular noise to behavioural variability in single bacterium Nature 428 574 
CrossRefPubMed

[13]

Lubashevsky I, Friedrich R and Heuer A 2009 Realization of levy walks as Markovian stochastic process Phys. Rev. E 79 011110 
CrossRef

[14]

Beck C 2001 Dynamical foundations of nonextensive statistical mechanics Phys. Rev. Lett. 87 180601 
CrossRef

[15]

Ausloos M and Lambiotte R 2006 Brownian particle having a fluctuating mass Phys. Rev. E 73 011105 
CrossRef

[16]

Alberts B, Johnson A, Walter P, Lewis J, Raff M and Roberts K 2007 Molecular Biology of the Cell 5th edn  (London: Garland Science) 

[17]

Wolfram S 1999 The Mathematica Book 4th edn  (Cambridge: Cambridge University Press) 

[18]

Luczka J, Talkner P and Hanggi P 2000 Diffusion of Brownian particles goverened by fluctuating friction Physica A 278 18–31 
CrossRef

[19]

Sher H and Montroll E W 1975 Anomalous transit-time dispersion in amorphous solids Phys. Rev. B 12 2455 
CrossRef

[20]

Mandelbrot B B and Van Ness J W 1968 Fractional brownian motions, fractional noises and applications SIAM Rev. 10 422–37 
CrossRef

[21]

Min W, Luo G, Cherayil B J, Kou S C and Xie X S 2005 Observation of a power-law memory kernel for fluctuations within a single protein molecule Phys. Rev. Lett. 94 198302 
CrossRefPubMed

[22]

Goychuk I and Hänggi P 2007 Anomalous escape governed by thermal 1/f noise Phys. Rev. Lett. 99 200601 
CrossRefPubMed

[23]

Burov S and Barkai E 2008 Critical exponent of the fractional Langevin equation Phys. Rev. Lett. 100 070601 
CrossRefPubMed

[24]

Muralidhar R, Ramkrishna D, Nakanishi H and Jacobs D 1990 Anomalous diffusion: a dynamic perspective Physica A 167 539–59 
CrossRef

[25]

Bouchaud J P and Georges A 1990 Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications Phys. Rep. 195 127 
CrossRef

[26]

Bel G and Barkai E 2005 Ergodicity breaking in continuous time random walk Phys. Rev. Lett. 94 240602 
CrossRef

[27]

Tolic-Norrelykke I M et al 2004 Anomalous diffusion in living yeast cells Phys. Rev. Lett. 93 078102 
CrossRefPubMed

[28]

He Y, Burov S, Metzler R and Barkai E 2008 Random time-scale invariant diffusion and transport coefficients Phys. Rev. Lett. 101 058101 
CrossRefPubMed

[29]

Sokolov I M, Lubelski A and Klafter Y 2008 Nonergodicity mimics inhomogeneity in single particle tracking Phys. Rev. Lett. 100 250602 
CrossRefPubMed

[30]

Deng W and Barkai E 2009 Ergodic properties of fractional Brownian–Langevin motion Phys. Rev. E 79 011112 
CrossRef

[31]

Dreisigmeyer D W, Stajic J, Nemenman I, Hlavacek W S and Wall M E 2008 Determinants of bistability in induction of the Escherichia coli lac operon IET Syst. Biol. 2 293–303 
CrossRefPubMed