Time-averaged MSD for switching diffusion

TL;DR

This paper analyzes a two-state switching diffusion model, deriving exact expressions for the mean and variance of the TAMSD, revealing ergodicity breaking effects in finite measurement times that align with experimental observations.

Contribution

It provides the first exact calculations of TAMSD statistics for switching diffusions, highlighting conditions for ergodicity breaking and offering an ergodic alternative to non-ergodic models.

Findings

Exact mean and variance of TAMSD derived.

Ergodicity breaking parameter can be arbitrarily large.

Switching diffusions can mimic non-ergodic behavior.

Abstract

We consider a classic two-state switching diffusion model from a single-particle tracking perspective. The mean and the variance of the time-averaged mean square displacement (TAMSD) are computed exactly. When the measurement time (i.e., the trajectory duration) is comparable to or smaller than the mean residence times in each state, the ergodicity breaking parameter is shown to take arbitrarily large values, suggesting an apparent weak ergodicity breaking for this ergodic model. In this regime, individual random trajectories are not representative while the related TAMSD curves exhibit a broad spread, in agreement with experimental observations in living cells and complex fluids. Switching diffusions can thus present, in some cases, an ergodic alternative to commonly used and inherently non-ergodic continuous-time random walks that capture similar features.

Equations

$$\overline{{\delta }^2}(t,T)=\frac{1}{T-t}\underset{0}{\overset{T-t}{\int }}d{t}_0(X(t+t_0)-X(t_0){)}^2,$$

$$\chi =\frac{⟨[\overline{{\delta }^2}(t,T){]}^2⟩}{⟨\overline{{\delta }^2}(t,T){⟩}^2}-1=\frac{\mathrm{var}\{\overline{{\delta }^2}(t,T)\}}{(\mathrm{mean}\{\overline{{\delta }^2}(t,T)\}{)}^2},$$

$$\frac{\partial }{\partial t}{P}_{1,i_0}$$

$$\frac{\partial }{\partial t}{P}_{2,i_0}$$

$$P_{i{i}_0}(x,t|x_0,t_0)=\underset{-\infty }{\overset{\infty }{\int }}\frac{dq}{2\pi }{e}^{-iq(x-x_0)}{\hat{P}}_{i{i}_0}(q,t-t_0),$$

$${\bf P}(q,t)=\left(\begin{array}[]{c c}\hat{P}_{11}(q,t)&\hat{P}_{12}(q,t)\\ \hat{P}_{21}(q,t)&\hat{P}_{22}(q,t)\\ \end{array}\right)=\exp(-{\bf M}_{q}t),$$

$${\bf M}_{q}=\left(\begin{array}[]{c c}D_{1}q^{2}+k_{12}&-k_{21}\\ -k_{12}&D_{2}q^{2}+k_{21}\\ \end{array}\right).$$

$${\bf M}_{q}{\bf V}_{q}={\bf V}_{q}\left(\begin{array}[]{c c}\gamma_{q}^{+}&0\\ 0&\gamma_{q}^{-}\\ \end{array}\right).$$

$${\gamma }_q^{\pm }$$

$${\bf V}_{q}=\left(\begin{array}[]{c c}\displaystyle\frac{k_{21}}{k_{12}+D_{1}q^{2}-\gamma^{+}_{q}}&\displaystyle\frac{k_{21}}{k_{12}+D_{1}q^{2}-\gamma^{-}_{q}}\\ 1&1\\ \end{array}\right).$$

$${\bf P}(q,t)=\frac{\left(\begin{array}[]{c c}e^{-\gamma_{q}^{+}t}\mu_{q}^{-}-e^{-\gamma_{q}^{-}t}\mu_{q}^{+}&k_{21}(e^{-\gamma_{q}^{-}t}-e^{-\gamma_{q}^{+}t})\\ k_{12}(e^{-\gamma_{q}^{-}t}-e^{-\gamma_{q}^{+}t})&e^{-\gamma_{q}^{-}t}\mu_{q}^{-}-e^{-\gamma_{q}^{+}t}\mu_{q}^{+}\\ \end{array}\right)}{\gamma_{q}^{+}-\gamma_{q}^{-}},$$

$$P(x_n,i_n,t_n;x_{n-1},i_{n-1},t_{n-1};\dots ;x_1,t_1,i_1;x_0,i_0,0)$$

$$=P_{i_n,i_{n-1}}(x_n,t_n|x_{n-1},t_{n-1})\dots P_{i_1,i_0}(x_1,t_1|x_0,0)$$

$$=\underset{{\mathbb{R}}^n}{\int }\frac{d{q}_1}{2\pi }\dots \frac{d{q}_n}{2\pi }{e}^{-i{q}_n(x_n-x_{n-1})-\dots -i{q}_1(x_1-x_0)}$$

$$\times {\hat{P}}_{i_n,i_{n-1}}(q_n,t_n-t_{n-1})\dots {\hat{P}}_{i_1,i_0}(q_1,t_1).$$

$$P(x_n,t_n;x_{n-1},t_{n-1};\dots ;x_1,t_1;x_0,0)$$

$$=\underset{{\mathbb{R}}^n}{\int }\frac{d{q}_1}{2\pi }\dots \frac{d{q}_n}{2\pi }{e}^{-i{q}_n(x_n-x_{n-1})-\dots -i{q}_1(x_1-x_0)}$$

$$\times {\mathcal{P}}_n(q_n,t_n-t_{n-1};\dots ;q_2,t_2-t_1;q_1,t_1),$$

$${\mathcal{P}}_n(q_n,t_n-t_{n-1};\dots ;q_2,t_2-t_1;q_1,t_1)$$

$$\displaystyle=\left(\begin{array}[]{c}1\\ 1\\ \end{array}\right)^{\dagger}{\bf P}(q_{n},t_{n}-t_{n-1})\ldots{\bf P}(q_{1},t_{1})\left(\begin{array}[]{c}p_{1}\\ p_{2}\\ \end{array}\right).$$

$$G(q)$$

$$⟨(X(t+t_0)-X(t_0){)}^2⟩=-\underset{q\to 0}{\lim }\frac{{\partial }^{2}G(q)}{\partial q^2}$$

$$=2\bar{D}t+2\frac{(p_1{k}_{12}-p_2{k}_{21})(D_1-D_2)}{k^2}(1-e^{-kt})e^{-k{t}_0},$$

$$k=k_{12}+k_{21}$$

$$\bar{D}=\frac{D_1{k}_{21}+D_2{k}_{12}}{k}$$

$$⟨X^2(t)⟩=2\bar{D}t+2\frac{(p_1{k}_{12}-p_2{k}_{21})(D_1-D_2)}{k^2}(1-e^{-kt}).$$

$$⟨\overline{{\delta }^2}(t,T)⟩=2\bar{D}t$$

$$+\frac{2(p_1{k}_{12}-p_2{k}_{21})(D_1-D_2)}{k^3(T-t)}(1-e^{-kt})(1-e^{-k(T-t)}).$$

$$p_1=p_1^{\mathrm{eq}}=\frac{k_{21}}{k},p_2=p_2^{\mathrm{eq}}=\frac{k_{12}}{k},$$

$$⟨\overline{{\delta }^2}(t,T)⟩=2\bar{D}t.$$

$$⟨[\overline{{\delta }^2}{]}^2⟩$$

$$\mathrm{var}\{\overline{{\delta }^2}\}=⟨[\overline{{\delta }^2}{]}^2⟩-⟨\overline{{\delta }^2}{⟩}^2.$$

$$G(q,q^{\prime })\equiv ⟨e^{iq(X(t_0+t)-X(t_0))+i{q}^{\prime }(X(t_0^{\prime }+t)-X(t_0^{\prime }))}⟩.$$

$$G(q,q^{\prime })={\mathcal{P}}_4(q^{\prime },t;0,t_0^{\prime }-t_0-t;q,t;0,t_0);$$

$$G(q,q^{\prime })={\mathcal{P}}_4(q^{\prime },t_0^{\prime }-t_0;q+q^{\prime },t_0+t-t_0^{\prime };q,t_0^{\prime }-t_0;0,t_0).$$

Full text

switching diffusion, Kärger model, ergodicity breaking, time averaged MSD

I Introduction

A continuous-time random walk (CTRW) was originally proposed to model hopping processes in semiconductors with randomly distributed heavy tailed waiting times due to random energetic trapping Montroll65 ; Bouchaud90 . This model has been extensively studied and rapidly became an archetypical model of anomalous diffusion in different contexts, including the intracellular transport in microbiology Metzler00 ; Bressloff13 ; Metzler14 . This is also an emblematic model of weak ergodicity breaking (WEB) and aging phenomena Bouchaud92 ; Burov11 ; Schulz13 . In fact, when the waiting time distribution has no first moment, there is no characteristic time scale for waiting events, so that the duration of a single stalling period may be comparable to the overall measurement time that prohibits self-averaging needed for ergodic behavior. Moreover, longer measurement times favor longer waiting times so that the statistical properties of the system depend on the measurement time. In the biological context, CTRWs were applied to mimic molecular caging in an overcrowded intracellular environment when a diffusing macromolecule, surrounded by other macromolecules and filaments, may wait a relatively long time in such an effective molecular “cage” before being able to jump to a next cage Wong04 ; Weigel11 ; Jeon11 ; Barkai12 . However, the validity of the considered heavy-tailed distribution of waiting times in the biological context remains debatable Szymanski09 . In particular, such a distribution may have an exponential cut-off so that anomalous features are transient, before getting into normal diffusion at long times Burov11 ; Lanoiselee16 . Although a CTRW with cut-off sounds more realistic, the underlying mathematical formalism is rather difficult.

In this light, a simpler model of random switching between several diffusion states, characterized by diffusion coefficient $D_{1},\ldots,D_{J}$, is appealing. A particle started in a state $i$ undergoes normal diffusion with diffusivity $D_{i}$ for a random time, until it switches to another state $j$, with the switching rate $k_{ij}$, and so on. Diffusion states can represent different conformations of a large macromolecule (e.g., globular versus filamentous structures) and thus distinct effective radii. Alternatively, diffusion states can account for eventual temporal binding of the diffusing particle to other molecules that may either slow down or even stop its motion. Such switching diffusions are often employed to describe the dynamics in biological systems Bressloff17 ; Sungkaworn17 ; Weron17 and used as simple models of intermittent processes Lanoiselee17 .

In this paper, we mainly focus on the two-state model with diffusion coefficients $D_{1}$ and $D_{2}$ and the switching rates $k_{12}$ and $k_{21}$. Note that $1/k_{12}$ and $1/k_{21}$ are the mean residence times of the particle in the states $1$ and $2$, respectively. The molecular caging effect can be modeled by setting $D_{2}=0$, i.e., the particle does not move in the second state. This model has been extensively studied, in particular, in the nuclear magnetic resonance literature, in which it is known as the Kärger model Karger85 ; Fieremans10 ; Moutal18 . Here, we look at this model from a single-particle tracking perspective. This simple model will allow us to investigate the reproducibility of measurements over individual random trajectories and the effect of their duration $T$. In particular, if both residence times $1/k_{12}$ and $1/k_{21}$ are small as compared to $T$, the particle switches very often between two states and manages to probe these states reliably during the measurement time. In this limit, the intermittent process is seen as an ergodic normal diffusion with some mean diffusivity $\bar{D}$ (see below). In contrast, if both residence times are much larger than $T$, the particle remains within a single state over the whole measurement with a high probability. As a consequence, such a single trajectory would bring information only about one state, whereas another trajectory may bring information only about the other state. In other words, individual trajectories are not representative of the whole dynamics and thus, from a practical point of view, not ergodic. Since the switching model is ergodic, we call this regime apparent WEB due to insufficient measurement time. Finally, the situation when one or both residence times are comparable to the measurement time is a borderline case. This model will thus help us to study the transition from one regime to the other that is quite typical for intracellular transport measurements.

In order to quantify the reliability of measurements and an eventual apparent WEB, we investigate the time-averaged mean square displacement (TAMSD) of a particle undergoing such an intermittent motion:

[TABLE]

where $t$ is the lag time. We compute exactly both the mean and the variance of the random variable $\overline{\delta^{2}}(t,T)$ to evaluate the ergodicity breaking (EB) parameter $\chi$, also known as the squared coefficient of variation of the TAMSD He08 ; Burov10 :

[TABLE]

where $\langle\ldots\rangle$ denotes the expectation, i.e., the ensemble average over the space of all possible trajectories $X(t)$ of duration $T$.

The paper is organized as follows. In Sec. II, we provide the propagator for multiple successive points in a matrix form which is then used to derive the mean and the variance of the TAMSD. Different asymptotic regimes of the resulting EB parameter are then analyzed. In Sec. III, we discuss an apparent WEB when the measurement time is comparable to or smaller than the mean residence times. In particular, we illustrate a large spread of TAMSD curves obtained via numerical simulations. We also mention extensions for multi-state switching diffusion and restricted diffusion in bounded domains. Section IV summarizes the results and concludes.

II Theoretical results

II.1 Propagator

We consider a particle that diffuses on a real axis ${\mathbb{R}}$ and switches randomly between two diffusivities $D_{1}$ and $D_{2}$ with rates $k_{12}$ and $k_{21}$ (see Sec. III.2 and III.3 for extensions). Let $P_{ii_{0}}(x,t|x_{0},t_{0})$ be the propagator of this particle, i.e., the probability density of finding the particle at point $x$ in state $i$ at time $t$ given that it started from point $x_{0}$ in state $i_{0}$ at time $t_{0}$. These four propagators satisfy four coupled partial differential equations

[TABLE]

(with $i_{0}=1,2$), subject to the initial conditions: $P_{i,i_{0}}(x,t=t_{0}|x_{0},t_{0})=\delta_{i,i_{0}}\delta(x-x_{0})$ (a rigorous mathematical formulation of switching models and of some their properties can be found in Yin10 ; Yin10b ; Baran13 ).

The solution of these equations is well known (e.g., see Karger85 ). Here we recall the main formulas that will be needed for the analysis of the TAMSD. In the Fourier space,

[TABLE]

Eqs. (3) are reduced to coupled first-order differential equations whose solution is easily found in a matrix form

[TABLE]

where

[TABLE]

The eigenvalues and eigenvectors of the matrix ${\bf M}_{q}$,

[TABLE]

are known explicitly:

[TABLE]

and

[TABLE]

One gets thus

[TABLE]

with $\mu_{q}^{\pm}=D_{1}q^{2}+k_{12}-\gamma_{q}^{\pm}$. In the special case $q=0$, one has $\gamma_{0}^{+}=k_{12}+k_{21}$ and $\gamma_{0}^{-}=0$.

To compute the moments of the TAMSD, one needs the propagator for multiple successive points, which is simply the product of the above propagators due to the Markov property of the process:

[TABLE]

Denoting by $p_{1}$ (resp., $p_{2}=1-p_{1}$) the probability of starting in the state $1$ (resp. 2) at time [math], the marginal propagator averaged over the state variables $i_{k}$ reads

[TABLE]

where

[TABLE]

II.2 Mean TAMSD

We first calculate the characteristic function of the displacement between times $t_{1}$ and $t_{2}$ such that $0<t_{1}<t_{2}$:

[TABLE]

From this characteristic function, one can compute the moments of the displacement, in particular, the mean square displacement:

[TABLE]

where

[TABLE]

and

[TABLE]

is the mean diffusivity. In particular, the ensemble averaged MSD is

[TABLE]

From Eq. (19), one can also deduce the mean TAMSD

[TABLE]

If the initial probabilities $p_{1}$ and $p_{2}$ are set to be from the equilibrium,

[TABLE]

then

[TABLE]

One cannot therefore reveal the intermittent character of this process from the mean TAMSD alone. Moreover, the mean value does not depend on the measurement time $T$, as for normal Brownian motion.

II.3 Variance of TAMSD

The computation of the variance of the TAMSD is much more involved as it requires the four-point correlation function. In fact, one has

[TABLE]

from which the variance follows as usual:

[TABLE]

For computing the second moment of the TAMSD, we introduce the characteristic function

[TABLE]

Since we have $t_{0}<t^{\prime}_{0}$ in Eq. (II.3), there are two cases:

(i) for $t_{0}<t_{0}+t<t^{\prime}_{0}<t^{\prime}_{0}+t$, we get

[TABLE]

(ii) for $t_{0}<t^{\prime}_{0}<t_{0}+t<t^{\prime}_{0}+t$, we get

[TABLE]

In the evaluation of the integral in Eq. (II.3), we consider separately two cases: $t<T/2$ and $t>T/2$.

Case $t<T/2$

In this case, one can split the integral in Eq. (II.3) into three parts

[TABLE]

(note that the integration variables $t_{0}$ and $t^{\prime}_{0}$ were shifted), with

[TABLE]

After long and cumbersome computations of these derivatives and integrals in Eq. (31), we derive the following expression for the variance of the TAMSD for $t<T/2$:

[TABLE]

where the initial probabilities $p_{i}$ were set to their equilibrium values $p_{i}^{\rm eq}$ in Eq. (24) to get a more compact expression.

Case $t>T/2$

In this case, only the second option in Eq. (30) is possible so that one does not need to split the integral in Eq. (II.3) into three parts, and one gets

[TABLE]

where the integration variable was shifted in the second integral. The evaluation of this integral yields the variance of the TAMSD for $t>T/2$:

[TABLE]

where the initial probabilities $p_{i}$ were again set to their equilibrium values $p_{i}^{\rm eq}$ for simplicity.

II.4 Ergodicity breaking parameter

The expressions (II.3, II.3) for the variance of the TAMSD present the main computational result of this paper. The first term in these expressions is the variance of the TAMSD for Brownian motion with the mean diffusivity $\bar{D}$ Qian91 ; Grebenkov11 ; Andreanov12 ; Grebenkov13 . One can check that the second term in Eqs. (II.3, II.3) is nonnegative, i.e., the switching between two states can only increase the variance of the TAMSD.

As Eqs. (II.3, II.3) for the variance are provided exclusively for the case $p_{i}=p_{i}^{\rm eq}$, for which the mean TAMSD in Eq. (25) is particularly simple, the analysis of the variance is equivalent, up to a simple multiplicative factor $(2\bar{D}t)^{2}$, to that of the ergodicity breaking parameter $\chi$ defined by Eq. (2). In the rest of the paper, we focus on this parameter.

As the trajectory duration $T$ goes to infinity (for a fixed lag time $t$), the EB parameter vanishes asymptotically as

[TABLE]

so that the switching process is ergodic, as expected. In contrast, in the double limit $k_{21}=p_{1}^{\rm eq}k\to 0$ and $k_{12}=p_{2}^{\rm eq}k\to 0$ (with fixed $p_{i}^{\rm eq}$ and $k\to 0$), Eq. (II.3) yields

[TABLE]

In this limit, the particle stays infinitely long in either of two states, i.e., the process is not ergodic, and the variance of the TAMSD does not vanish in the limit $T\to\infty$. This singular situation can also describe two populations of particles with distinct diffusivities $D_{1}$ and $D_{2}$, and the leading term of the EB parameter in Eq. (36) is a consequence of their mixture with relative fractions $p_{1}^{\rm eq}$ and $p_{2}^{\rm eq}$. This limit could alternatively be obtained by setting $\overline{\delta^{2}}=\alpha\,\zeta_{1}+(1-\alpha)\,\zeta_{2}$, where $\zeta_{i}$ is the TAMSD for the $i$-th population (that differs by the factor $D_{i}$), and a random selection between two populations is realized by a Bernoulli random variable $\alpha$ taking the value $1$ with probability $p_{1}^{\rm eq}$ and [math] with probability $p_{2}^{\rm eq}=1-p_{1}^{\rm eq}$. Moreover, the EB parameter $\chi$ remains close to the limiting expression (36) when $T\ll 1/k$. In other words, if the measurement time $T$ is short as compared to the residence time $1/k$, the system exhibits an apparent WEB. Clearly, the order of two limits, $k_{12},\,k_{21}\to 0$ and $T\to\infty$, does matter here: sending $T\to\infty$ for fixed $k_{12}$ and $k_{21}$ yields the zero variance, as expected.

In the limit $t\to 0$, Eq. (II.3) gives

[TABLE]

We first consider the particular case, in which two switching rates are identical: $k_{12}=k_{21}=k/2$. For our illustrative purposes, we use dimensionless units for all parameters. We set the measurement time $T=1000$ (i.e., the trajectory with a thousand steps). We recall that the EB parameter for Brownian motion is a monotonously growing function of the lag time $t$, so that the smallest available lag time $t=1$ provides the most accurate estimation of the TAMSD (see, e.g., Grebenkov11 ). One can check that the same property holds for two-state switching diffusion. As a consequence, we select the lag time $t=1$ to be in the optimal situation. Figure 1 shows the behavior of the EB parameter as a function of the residence time $1/k$. In the case of equal diffusion coefficients, $D_{2}/D_{1}=1$, the states are identical, and the switching model is reduced to normal diffusion. The EB parameter does not depend on the switching rate and is equal (in the leading order) to $4t/(3T)$ according to Eq. (II.4). In turn, the stronger the difference between two states (i.e., smaller $D_{2}/D_{1}$), the larger the EB parameter at large residence times $1/k$. In the extreme case of one immobile state (with $D_{2}=0$), the EB parameter reaches the value $1$ according to Eq. (36). In contrast, a fast switching even between very distinct states leads again to normal diffusion with the mean diffusivity $\bar{D}$ and thus the EB parameter remains close to $4t/(3T)$. Note that in the considered case of equal switching rates, the EB parameter does not exceed $1$. We conclude that, due to slow switching between states, the distribution of TAMSD even for the smallest lag time ($t=1$) can be relatively broad, i.e., the standard deviation can be comparable to the mean value: $\chi\sim 1$. In this situation, an increase of the trajectory duration $T$ can improve the estimation only when $T$ exceeds the mean residence time $1/k$. This is in contrast with the case of normal diffusion, for which the EB parameter is of the order of $t/T$ and can thus be made very small, allowing for accurate estimations of the diffusion coefficient for long enough trajectories.

A richer insight on the EB parameter is provided in Fig. 2 that shows the dependence on two residence times $1/k_{12}$ and $1/k_{21}$. In this situation, there are four relevant time scales: $t$, $T$, $1/k_{12}$ and $1/k_{21}$ and thus various regimes. As previously, we fix $t=1$ and $T=1000$. First, we present in Fig. 2(a) the case of an immobile state with $D_{2}=0$ (and $D_{1}=1$). When $1/k_{21}$ is the smallest time scale, the particle does not almost stay in the immobile state, and the EB parameter is close to its value $4t/(3T)$ for normal diffusion. When $1/k_{21}$ is getting comparable to the lag time $t$, the switching is not fast enough any more, and the presence of the immobile phase increases the EB parameter. Finally, when $1/k_{21}$ increases further, the EB parameter can become much larger than $1$. In fact, in the limit $k_{21}\to 0$ (with fixed $k_{12}$), one gets

[TABLE]

where

[TABLE]

One can see that the EB parameter can be made arbitrarily large by decreasing $k_{21}$. Indeed, a particle started in the immobile state mainly remains in this state, whereas a particle started in the mobile state becomes immobile with the switching rate $k_{12}$. As a consequence, random trajectories may have a broad distribution of stalling periods and thus the broad distribution of TAMSD. We note that the other mean residence time, $1/k_{12}$, also influences the EB parameter but it is less relevant than $1/k_{21}$.

The situation is considerably different for a particle undergoing slow diffusion in the second state (i.e., $D_{2}$ is small but not strictly zero). Figure 2(b) shows an example for $D_{2}=0.01$. While the behavior of the EB parameter for small $1/k_{21}$ is expectedly similar to the former case with $D_{2}=0$, there is significant difference for large $1/k_{21}$. First, the values of the EB parameter, which can still be large, are much smaller than those shown in Fig. 2(a). Second, the EB parameter exhibits a maximum as a function of the mean residence time $1/k_{12}$ for a fixed $1/k_{21}$. When $1/k_{12}$ is small (with $1/k_{21}$ large), the particle stays most of the time in the second state and thus undergoes almost normal diffusion with diffusivity $D_{2}$, so that the EB parameter again reaches the small value $4t/(3T)$. This is a significant difference with respect to the case $D_{2}=0$.

III Discussion

When the measurement time $T$ is comparable to or smaller than the residence times $1/k_{12}$ and $1/k_{21}$, the particle does not have enough time to probe the phase space. An individual random trajectory of the particle is thus not representative of the ensemble, suggesting an apparent weak ergodicity breaking. In particular, the EB parameter can be of order of $1$ and even much larger, particularly if the diffusion coefficient in one of the states is very small (or zero). As a consequence, one can expect a large spread of TAMSD curves evaluated from individual trajectories. This feature, which is often observed in experiments (see, e.g., Golding06 ; Jeon11 ; Mosqueira18 ), was often interpreted as an indication of non-ergodicity and thus attributed to CTRW as a basic non-ergodic model. Here we showed that the ergodic two-state model can lead to similar features.

III.1 Spread of TAMSD curves

In order to illustrate the spread of TAMSD curves, we perform Monte Carlo simulations. For a prescribed set of parameters $D_{1}$, $D_{2}$, $k_{12}$ and $k_{21}$, we set the trajectory duration $T=1000$ (i.e., the number of steps) and generate $T$ random centered Gaussian increments with unit variance. We also generate a sequence of successive residence times in two alternating states according to the exponential laws with the rates $k_{12}$ and $k_{21}$, that results in a random sequence of state variables $i_{1},i_{2},\ldots$ (each $i_{n}$ taking values $1$ or $2$). All the increments are rescaled by $\sqrt{2D_{i_{n}}}$ and then cumulatively summed up to produce a random trajectory, from which the TAMSD is computed for all lag times from $1$ to $T-1$ by discretizing Eq. (1). This computation is repeated $10\,000$ times to obtain a reliable statistics. Fixing $D_{1}=1$, we are left with $D_{2}$, $k_{12}$ and $k_{21}$ as the major parameters.

Figure 3 illustrates the spread of TAMSD curves for three sets of parameters. The first set with $D_{2}=1$ corresponds to normal diffusion (here, switching rates do not matter as $D_{1}=D_{2}$). In this case, TAMSD curves are close to each other at small lag times and then getting more spread at larger lag times because the time average becomes less and less efficient (Fig. 3a). In the second set (with $D_{2}=0.01$, $k_{12}=k_{21}=0.1$), the mean residence times $1/k_{12}$ and $1/k_{21}$ are chosen to be much smaller than the measurement time $T$. As switching is rapid enough, TAMSD curves remain close to each (Fig. 3b), as for normal diffusion. In the third set of parameters, we keep $D_{2}=0.01$ but decrease both switching rates: $k_{12}=10^{-2}$ and $k_{21}=10^{-3}$. While $1/k_{12}$ is still much smaller than $T$, $1/k_{21}$ is equal to $T$ that leads to a wide spread of TAMSD curves (Fig. 3c), in agreement with large values of the EB parameter in this case. Similar spreads were observed in single-particle tracking experiments in living cells (see, e.g., Golding06 ; Jeon11 ; Mosqueira18 ).

Another way of presenting the spread consists in plotting the distribution of TAMSD for a fixed lag time. The distribution of TAMSD for ergodic Brownian motion and other Gaussian processes was studied in Grebenkov11 ; Andreanov12 ; Sikora17 ; Gajda18 , while the analysis of its asymptotic behavior for non-ergodic CTRW was initiated in He08 ; Lubelski08 (see also reviews Burov11 ; Metzler14 ). As our computation for two-state switching diffusion is limited to the first two moments, we show in Fig. 4 the empirical distribution of TAMSD obtained from simulated trajectories for the same three sets of parameters. The empirical distributions are presented for three lag times: $t=1$, $t=10$ and $t=105$. As expected, the distribution is getting larger with the lag time, reflecting less and less efficient time averaging. For the second set of parameters, the distributions are close to that for normal diffusion (compare Figs. 4a and 4b). In contrast, the distribution for the third set of parameters is much broader and almost does not depend on the lag time. This is the characteristic feature of non-ergodic dynamics (e.g., the distribution of TAMSD for CTRW is also broad).

To complete this discussion, we outline a subtle difference in the behavior of the EB parameter between a continuous-time process and its discrete-time approximation. Let us first discuss Brownian motion, which is often approximated by random walks on a lattice or by a sequence of Gaussian jumps. While such approximations are known to converge to Brownian motion (see, e.g., Morters ; Grebenkov16 ), some functionals involving integrals along sample trajectories of these processes can be different. So, the variance of the TAMSD for continuous-time Brownian motion, given by the first term in Eqs. (II.3, II.3), is different from the variance of the TAMSD for a discrete-time random walk with Gaussian jumps, first derived by Qian et al. Qian91 and later analyzed in Grebenkov11 ; Grebenkov13 . In particular, when $t\ll T$, the leading term of the EB parameter is $4t/(3T)$ in the continuous case, and $2t/T$ in the discrete case. While the scaling form, $t/T$, is identical for both cases, the numerical prefactor is different. This difference is related to the fact that an erratic trajectory of Brownian motion between two successive discrete times is replaced by a jump in the discrete case.

There is a similar distinction between the continuous-time switching diffusion, studied in Sec. II, and discrete-time Monte Carlo simulations presented here. In particular, the variance of the TAMSD, computed via the exact formulas (II.3, II.3), would differ from its numerical estimation by Monte Carlo simulations. In spite of this subtle difference between continuous-time and discrete-time processes, the qualitative conclusion on distinct diffusivity states and insufficient measurement time as eventual causes of the apparent weak ergodicity breaking remains valid, as confirmed by numerical results of this section. The exact computation of the EB parameter for a discrete-time switching diffusion by adapting the general technique from Grebenkov11 presents an interesting perspective for future research.

III.2 Extension to multi-state models

The above computation can be formally extended to a switching model with $J$ states, which is characterized by a set of diffusion coefficients $D_{i}$ and switching rates $k_{ij}$. The propagators satisfy for each $i_{0}=1,\ldots,J$:

[TABLE]

where $k_{ii}\equiv-\sum\limits_{j\neq i}k_{ij}$. The formula (12) for the marginal propagator for multiple successive points in the Fourier space remains practically unchanged:

[TABLE]

where $p_{i}$ is the probability of starting in the state $i$, ${\bf P}(q,t)=\exp(-{\bf M}_{q}t)$, and ${\bf M}_{q}$ is now the $J\times J$ matrix of the form ${\bf M}_{q}=q^{2}{\bf D}-{\bf K}^{\dagger}$, with $({\bf D})_{ij}=\delta_{ij}D_{i}$ and $({\bf K})_{ij}=k_{ij}$. As the eigenvalues $\gamma_{q}$ of the matrix ${\bf M}_{q}$ are determined as the zeros of the polynomial

[TABLE]

there is no explicit formula in general. However, as we are interested in evaluating the limit $q\to 0$ to get the mean square displacement and other related quantities, one can apply the standard perturbation theory by treating $q^{2}{\bf D}$ as a perturbation of the matrix ${\bf M}_{0}=-{\bf K}^{\dagger}$. As the matrix ${\bf M}_{0}$ is neither symmetric, nor invertible, the analysis is more involved. While analytical calculations for a general multi-state model seem challenging, numerical computations of the propagator and related quantities (such as the mean square displacement) are efficient due to the matrix form (40) when the number of states is not too high (say, below a thousand). It is also worth noting that multi-state switching diffusions can be seen as a discrete version of diffusing diffusivity models Chubynsky14 ; Jain16 ; Chechkin17 ; Lanoiselee18a ; Lanoiselee18b ; Sposini18 ; Sposini19 , in which the diffusivity $D_{t}$ of a tracer changes continuously in time (see Grebenkov19 for details).

III.3 Extensions to switching diffusion in bounded domains

While we focused on one-dimensional diffusion, the above computation can be carried on for more sophisticated processes. Here, we briefly mention switching diffusion in an arbitrary bounded domain $\Omega\subset{\mathbb{R}}^{d}$, for which the propagator for multiple successive points, $P(\bm{x}_{n},i_{n},t_{n};\,\bm{x}_{n-1},i_{n-1},t_{n-1};\,\ldots;\bm{x}_{0},i_{0},0)$, and its average over the state variables $i_{k}$, $P(\bm{x}_{n},t_{n};\,\bm{x}_{n-1},t_{n-1};\,\ldots;\bm{x}_{0},0)$, can be expressed in analogy to Eqs. (10, 11). For instance, Eq. (11) becomes

[TABLE]

where ${\mathcal{P}}_{n}$ is given by Eq. (12) (or Eq. (40) in a multi-state case), while $\lambda_{n}$ and $u_{k}(\bm{x})$ are the eigenvalues and $L_{2}$-normalized eigenfunctions of the Laplace operator $\Delta$, satisfying $\Delta u_{k}(\bm{x})+\lambda_{k}u_{k}(\bm{x})=0$ in $\Omega$. Boundary conditions determine the properties of the boundary in a standard way Redner : Neumann condition incorporates passive impermeable walls whereas Dirichlet or Robin conditions allow one to describe diffusion-controlled reactions and the related first-passage time statistics, see Lanoiselee18b (note that ${\mathcal{P}}_{1}(q,t)$ corresponds to $\Upsilon(t;q^{2})$ in Lanoiselee18b ). For instance, one can investigate diffusion-controlled reactions in presence of buffers that can reversely bind the diffusing molecule and thus affect its mobility. We emphasize that the above derivation requires that switching modifies only the diffusivity of the particle but does not change other properties (e.g., reactivity). In this case, restricted diffusion in each state is characterized by the same set of Laplacian eigenmodes. In contrast, this extension is not applicable to other intermittent processes Benichou11b such as, e.g., surface-mediated diffusion Benichou10 ; Benichou11 ; Rojo11 ; Rupprecht12a ; Rupprecht12b , two-channel diffusion with switching reactivity Godec17 or switching boundary conditions Bressloff15 ; Lawley16 .

From the propagator, one can compute the mean square displacement, as well as the mean and the variance of the TAMSD. However, these formulas include multiple infinite sums over eigenvalues and thus remain formal, so we do not present these results.

IV Conclusion

We revisited the two-state switching diffusion model from a single-particle tracking perspective. The exact formulas for the mean and the variance of the TAMSD were derived, in order to investigate the behavior of the ergodicity breaking parameter. When the mean residence times $1/k_{12}$ and $1/k_{21}$ are small as compared to the measurement time $T$, switching is fast enough for a particle to probe both states so that a single trajectory is representative of the ensemble. In particular, the distribution of the TAMSD is narrow to allow for an accurate estimation of the mean diffusivity $\bar{D}$. Roughly speaking, the two-state switching diffusion looks as normal diffusion with diffusivity $\bar{D}$. In contrast, when $1/k_{12}$ and/or $1/k_{21}$ are comparable to or larger than $T$, an individual trajectory is not representative of the ensemble, the EB parameter may exceed $1$, and the distribution of the TAMSD is broad. As a consequence, TAMSD curves exhibit a significant spread, in agreement with experimental observations in living cells. While the two-state switching diffusion is ergodic (when $k_{12}$ and $k_{21}$ are strictly positive), the analysis of individual trajectories may indicate non-ergodic features. As this fictitious non-ergodicity is the mere consequence of insufficient measurement time $T$, we called it “apparent weak ergodicity breaking”. A similar behavior was observed for CTRW with an exponential cut-off: when the measurement time is much smaller than the cut-off time, the trajectory “looks” as that of usual non-ergodic CTRW but when the measurement time grows and exceeds the cut-off time, the non-ergodic features disappear Lanoiselee16 . The relative simplicity of switching diffusions presents their advantage as potential models for describing some intracellular and on-membrane transport Sungkaworn17 ; Weron17 . Moreover, these are natural models for describing the effect of buffers that can reversely bind the diffusing particle and thus affect its mobility. More generally, our results illustrate that a broad spread of TAMSD curves, often observed in experiments, is not necessarily a signature of weak ergodicity breaking. It may just be a consequence of the insufficient duration of acquired trajectories. Statistical tests of ergodicity breaking need to be systematically performed in such situations Magdziarz11 ; Lanoiselee16 .

Acknowledgements.

DG acknowledges the financial support by French National Research Agency (ANR Project ANR-13-JSV5-0006-01).