Steady diffusion in a drift field: a comparison of large deviation techniques and multiple-scale analysis

TL;DR

This paper compares large deviation techniques and multiple-scale analysis to derive effective advection-diffusion properties of particles with internal states in a force field, providing new methods and insights into their diffusion behavior.

Contribution

It introduces two novel methods—large deviation and multiple-scale analysis—for deriving effective diffusion and drift in particles with internal states, applicable to both discrete and continuous cases.

Findings

Effective drift velocity is proportional to the averaged mobility.

Effective diffusion has an equilibrium term and a force-dependent quadratic term.

Measuring diffusion components can determine kinetic rates from moments.

Abstract

A particle with internal unobserved states diffusing in a force field will generally display effective advection-diffusion. The drift velocity is proportional to the mobility averaged over the internal states, or effective mobility, while the effective diffusion has two terms. One is of the equilibrium type and satisfies an Einstein relation with the effective mobility while the other is quadratic in the applied force. In this contribution we present two new methods to obtain these results, on the one hand using large deviation techniques, and on the other by a multiple-scale analysis, and compare the two. We consider both systems with discrete internal states and continuous internal states. We show that the auxiliary equations in the multiple-scale analysis can also be derived in second-order perturbation theory in a large deviation theory of a generating function (discrete internal states) or generating functional (continuous internal states). We discuss that measuring the two components of the effective diffusion give a way to determine kinetic rates from only first and second moments of the displacement in steady state.

Figures

Equations

$$D^{eff}=\frac{T\left(k_{+}{\mu }_{-}+k_{-}{\mu }_{+}\right)}{k_{+}+k_{-}}+\frac{F^2{\left({\mu }_{+}-{\mu }_{-}\right)}^2{k}_{+}{k}_{-}}{(k_{+}+k_{-}{)}^3}$$

$$d{x}_t=F{\mu }_{a}dt+\sqrt{2T{\mu }_a}d{W}_t$$

$$\frac{d}{dt}\left(\begin{array}[]{c }P_{+}\\ P_{-}\end{array}\right)=\left(\begin{array}[]{c c }-k_{+}&k_{-}\\ k_{+}&-k_{-}\end{array}\right)\left(\begin{array}[]{c }P_{+}\\ P_{-}\end{array}\right)=\mathsf{K}\left(\begin{array}[]{c }P_{+}\\ P_{-}\end{array}\right)$$

$$w_{+}=\frac{k_{-}}{k_{-}+k_{+}}{w}_{-}=\frac{k_{+}}{k_{-}+k_{+}}.$$

$$\frac{\partial P_a(x)}{\partial t}$$

$$\{a_0\text{mbox}in(0={\tau }_0,{\tau }_1);\cdots ;a_n\text{mbox}in({\tau }_n,{\tau }_{n+1});\cdots ;$$

$$\cdots ;a_N\text{mbox}in({\tau }_N,{\tau }_{N+1}=t)\}$$

$$P(\Delta x_n|a_n;\Delta {\tau }_n)=\frac{1}{\sqrt{4\pi \Delta {\tau }_{n}T{\mu }_{a_n}}}{e}^{-\frac{{\left(\Delta x_n-F{\mu }_{a_n}\Delta {\tau }_n\right)}^2}{4\Delta {\tau }_{n}T{\mu }_{a_n}}}.$$

$$P(\Delta x|\{a\};t)=\frac{1}{\sqrt{4\pi T{\sum }_n\Delta {\tau }_n{\mu }_{a_n}}}{e}^{-\frac{{\left(\Delta x-F{\sum }_n{\mu }_{a_n}\Delta {\tau }_n\right)}^2}{4T{\sum }_n\Delta {\tau }_n{\mu }_{a_n}}}.$$

$$P(\Delta x|\{a\};t)=P(\Delta x|f;t)=\frac{e^{-\frac{{\left(\Delta x-⟨\Delta x|f⟩\right)}^2}{4tD(f)}}}{\sqrt{4\pi tD(f)}}$$

$$⟨\Delta x|f⟩=Ft\mu (f)$$

$$D(f)\equiv T\mu (f)=\frac{T}{tF}⟨\Delta x|f⟩$$

$$\mu (f)=f{\mu }_{+}+(1-f){\mu }_{-}$$

$$P(\Delta x)=\int P(\Delta x|f)P(f)df.$$

$$G_s^a(t)=\overline{e^{s{\tau }_{+}}\delta (a_t=a)}$$

$$\frac{d}{dt}\left(\begin{array}[]{c }G^{+}_{s}\\ G^{-}_{s}\end{array}\right)=\tilde{\mathsf{M}}_{s}\left(\begin{array}[]{c }G^{+}_{s}\\ G^{-}_{s}\end{array}\right)=\left(\begin{array}[]{c c }-k_{+}+s&k_{-}\\ k_{+}&-k_{-}\end{array}\right)\left(\begin{array}[]{c }G^{+}_{s}\\ G^{-}_{s}\end{array}\right)\,.$$

$${\lambda }_0(s)=\frac{1}{2}\left(-k_{-}-k_{+}+s+\sqrt{(k_{-}+k_{+}-s{)}^2+4k_{-}s}\right).$$

$$\overline{{\tau }_{+}}$$

$$\overline{{\left({\tau }_{+}-\overline{{\tau }_{+}}\right)}^2}$$

$$\overline{{\left({\tau }_{+}-\overline{{\tau }_{+}}\right)}^3}$$

$$\overline{{\left({\tau }_{+}-\overline{{\tau }_{+}}\right)}^4}$$

$$⟨\Delta x⟩$$

$$⟨(\Delta x-⟨\Delta x⟩{)}^2⟩=$$

$$⟨(\Delta x-⟨\Delta x⟩{)}^3⟩$$

$$\displaystyle\frac{\partial}{\partial t}\left(\begin{array}[]{c }P_{+}(x)\\ P_{-}(x)\end{array}\right)$$

$$=$$

$${\mathcal{L}}_a^{†}$$

$$\tilde{x}={\epsilon }^{0}x,X={\epsilon }^{1}x$$

$$\theta ={\epsilon }^{0}t,\vartheta ={\epsilon }^{1}t,\tau ={\epsilon }^{2}t,$$

$$P=P_a(\theta ,\vartheta ,\tau ,\tilde{x},X).$$

$$\frac{\partial }{\partial t}$$

$$\frac{\partial }{\partial x_i}$$

$$P=P^{(0)}+\epsilon P^{(1)}+{\epsilon }^2{P}^{(2)}+\dots ,$$

Full text

I Introduction

Measuring on microscopic and mesoscopic scales is the art of rendering visible what we cannot directly see. If the phenomenon of interest is faster than your time resolution you will only see a blur, but what happens on the faster scale may still be crucial. For instance, the fraction of the time a protein can be found in a given state is proportional to the Gibbs factor, which is generally the ratio of an on rate and an off rate, but the speed of catalyzed reactions depend on the kinetic rates directly english2006 . Similarly, the average occupancy of transcription factors to their binding sites on DNA determines the level of gene expression, but the time scales of regulation of transcription are set by the on rates and off rates separately Buchler2003 ; Alon . State-of-the-art methods to measure kinetic rates rely on precise imaging of the molecules motion and translate dynamic changes of their diffusion into transition rates (see persson2013 ; persson2013review ; hammar2014 for in vivo measurements and terazima2006 for an in vitro method). These methods have had great impact in determining the actual time scales and mechanisms of gene regulation and elucidating classical issues in quantitative molecular biology such as how a transcription factor finds its binding sites on DNA BergvonHippel85 ; BergvonHippel89 ; HalfordMarko2004 ; hammar2014 . They are however challenging to use, and other approaches, even if of more limited utility, may also be of value.

In this paper we study a non-equilibrium effect, the dependence of the effective diffusion on an external applied force. The effect can hardly be said to be new. In fact, to our best knowledge it was first demonstrated by K.J. Mysels in a paper published a little more than 60 years ago Mysels56 . This author observed that if particle or system transits between two internal states as $+\rightleftharpoons-$ with rates $k_{+}$ and $k_{-}$, and if the mobilities of the particle or system in the two states are $\mu_{+}$ and $\mu_{-}$, then in a force field $F$ the effective diffusion is

[TABLE]

In the first term $T$ is the temperature in units such that Boltzmann’s constant $k_{B}$ is equal to unity. The second term, in Mysels56 called the “electrodiffusion coefficient”, is the a priori surprising term, which cannot simply be postulated by analogy to fluctuation-dissipation phenomena close to equilibrium. Similar terms were discussed quite some time ago by Van den Broeck and co-workers and by Mackey and co-workers, see e.g. VanDenBroeck83 ; VanDenBroeck84 ; Mackey97 , and recently reconsidered in the context of Taylor dispersion taylor1953 ; vandenbroeck1990 ; kahlen2017 . The effective diffusion of a protein was recently analyzed in illien2016diffusion and applied to fluorescence correlation spectroscopy (FCS) measurements in illien2017exothermicity . These investigations concern non-equilibrium phenomena where the protein functions as an enzyme acting on a substrate but they do not feature any externally applied force.

The first goal of the present paper is to present (1) from a modern perspective. We will derive it both by means of the scaled cumulant generating function and by a multiple-scale analysis. The cumulant technique also gives access to higher moments of the displacement. We will also show that a similar behavior as (1), which we recently found for the case of a rotating Brownian particle with a continuous internal state (the particle’s orientation) aurell2016 can also be derived from a large deviation argument and second order time-independent perturbation theory. The second goal of the paper is to bring effects such as (1) anew to the attention of the single-molecule biophysics community, as a possible means to indirectly measure kinetic rates. Although established a long time ago these non-equilibrium steady-state relations appear to have been somewhat occluded, and the recent revival of interest has been motivated by other questions. The analysis of Taylor dispersion in kahlen2017 is for instance partly foreshadowed by an analysis of active Brownian motion (Janus particles) pietzonka2016extreme using cumulant generating functions, from which expressions like Eq. (1) can be derived.

The paper is organized as follows. In Section II we define the model and in Section III we give a heuristic analysis that establishes a term like the second one in (1), with a quadratic dependence on the applied force. In Section IV we derive this result systematically using large deviation theory varadhan2008large ; touchette2009large , and in Section V we obtain the results by means of multiple-scales analysis pavliotis08 . In Section VI we generalize to space-dependent switching rates and to several internal states, details of the latter analysis given in Appendix A. In Section VII we discuss the effect as a device to measure kinetic rates and in Section VIII we sum up and discuss our results.

Some technical material for the explicit solutions of section  V is given in Appendix B.

In two additional Appendices C and D we extend for completeness the large deviation analysis to the case of a continuum of internal states, of which one example is a Brownian particle undergoing translation and rotation. We show that an auxiliary equation which appeared in an earlier multiple-scale analysis aurell2016 ; marino2016b , can also be derived by second-order perturbation theory of generating function of the time spent in different orientations. We further show that the large deviation analysis can also be pushed to the third centered moment which we show to generically increase linearly in time, see Appendix D.

II Model

To illustrate our method we consider a particle diffusing in a fluid at temperature $T$ subject to an external constant force $F$ (Boltzmann constant set to 1). The particle can switch between two different states ($-$ and $+$) associated with different mobilities ($\mu_{-}$ and $\mu_{+}$) as sketched in Figure 1. The system can be described by a Langevin equation

[TABLE]

where the internal variable $a$ can be in either of two states denoted $-$ and $+$ and $dW_{t}$ is a Wiener increment. The internal state changes Markovianly with exponentially distributed waiting times and with transition rates $k_{+}$ (from state $+$ to $-$) and $k_{-}$ (for the reversed transition $-$ to $+$) that are independent of spatial coordinate $x$. For simplicity we consider constant rates in time. The probability of finding the internal variable in a given state then evolves according to the master equation:

[TABLE]

and the steady state probability reads:

[TABLE]

The joint probability density of being in position $x$ and state $a$ is determined by mixed diffusion and Master equation:

[TABLE]

Our aim is to describe the statistics of the particle displacement on time scales much longer than the typical relaxation of the discrete internal state.

III Heuristic analysis

The gist of the idea can be explained in simple terms as follows. At the microscale, the mobility of a body determines how its motion responds to an applied force and, together with temperature sets the magnitude of the thermal fluctuations, see (2) above. If the particle dynamically switches between states with different mobilities, how far the body moves under a force will depend on the time spent in each state. The average displacement is determined by the magnitude of the external force and by the effective mobility which is set by the mobility in each internal state weighted by the average fraction of time spent in a state.

On top of such averages are fluctuations in how far the body travels. Part of these fluctuations are thermal, due to the collisions with the fluid molecules in which the body is diffusing, and therefore proportional to the effective mobility (Einstein relation). However, given the presence of the external force, there are also fluctuations caused by the switching between states of different mobility. Indeed, a particle that has been more (less) often than usual in the state with the larger friction will lag behind (be ahead of) the average position, and this increases the spread around the average position in a way that is analogous to the phenomenon known as Taylor dispersion taylor1953 ; vandenbroeck1990 ; kahlen2017 . Such contribution is proportional to the variance of the time spent in each state and is quadratic in the applied force. In the rest of this paper we will discuss this issue more quantitatively, determine the prefactor of the square of the force that depends on the kinetic rates, and compare and contrast two ways of carrying out the analysis. To illustrate the problem we will start with the simple case in which the internal variable can take only two states and then generalize it.

IV Large deviation analysis

In this section we will use the scaled cumulant generating function to turn the above qualitative argument into a quantitative prediction. The long term statistics of the particle displacement can be computed by considering the explicit expression of its probability. The dynamics of the discrete internal variable is not influenced by the position of the particle. Its trajectories from an initial time [math] to a final one $t$ are specified by the sequence of states $a(\tau)$ with $0\leq\tau\leq t$. Let us number the interval between transitions as $n$ and the transition times as $\tau_{n}$. We then have that in an interval $n$, of duration $\tau_{n+1}-\tau_{n}$, the system resides in state $a_{n}$

[TABLE]

as depicted in Figure 1. Between transitions, the internal variable is fixed at its value $a_{n}$ and the particle undergoes a Langevin dynamics with a fixed value of the mobility $\mu_{a_{n}}$. The probability density of observing a displacement $\Delta x_{n}$ in the time interval $\Delta\tau_{n}=\tau_{n+1}-\tau_{n}$ in which $a=a_{n}$ is a simple Gaussian:

[TABLE]

The spatial increments $\Delta x_{n}$ depend on the realizations of the Wiener noise in the corresponding time intervals $\Delta\tau_{n}$, and are therefore independent random variables. The probability of overall displacement $\Delta x=\sum_{n}\Delta x_{n}$ is hence also Gaussian

[TABLE]

The only information about the dynamics of the discrete internal variable affecting the displacement probability is the total time spent in one of the two states, e. g., $\tau_{+}=\sum_{n}\tau_{n}\delta_{a_{n}}^{+}$ where $\delta_{a_{n}}^{+}$ is an indicator function that is $1$ when the discrete state $a_{n}=+$ and [math] otherwise. Introducing the the fraction of time spent in the $+$ state: $f$ and recalling that $\sum_{n}\tau_{n}=t$, expression (7) simplifies to

[TABLE]

where we have introduced the average particle displacement conditioned on spending a fraction $f$ of the total time $t$ in state $+$

[TABLE]

and

[TABLE]

with

[TABLE]

which is linear in $f$. In principle, it is then possible to compute the marginal probability of the body displacement by integrating the conditional expression (8) multiplied by $P(f)$, the probability of the fraction of time spent in the $+$ state, also known as the empiric (sample) average of state $+$:

[TABLE]

Despite the apparent simplicity of the problem, the explicit expression for $P(f)$ is rather complicated (see e.g. giddings55 ) and the integral cannot be analytically computed in closed form. However, as we shall show below, given the Gaussianity of $P(\Delta x|f)$, the moments of the displacement $\Delta x$ can be computed in a straight-forward manner starting from the moments of $f$. These moments of the fraction for time spent in the $+$ state further simplify in the long time limit, and hence lead to Eq. (1).

Scaled cumulant generating function of the fraction of time spent in a state.

The moment generating function of the time spent in state $+$: $\tau_{+}$ conditional on being in state $a$ at time $t$ is given by

[TABLE]

where the overline represents the average with respect of the distribution $P(f)$ for a process of duration $t$. The change in time of the generating function can be found by the methods used e.g. in Gardiner (Section 7.1) to analyze Birth-Death Master Equations and is given by

[TABLE]

For long times the moment generating function will be dominated by the leading eigenvalue of the matrix $\tilde{\mathsf{M}}_{s}$, which, for the case under consideration, is

[TABLE]

The scaled cumulant generating function is defined as (see e.g. touchette2009large ) $\lim_{t\to\infty}\frac{1}{t}\log G_{s}=\lambda_{0}(s)$ and can be used to evaluate the long-time expressions of the cumulants of the time spent in a state:

[TABLE]

Higher-order cumulants can also be found from higher derivatives of $\lambda$ in (15).

Long-time moments of the displacement.

The above explicit expressions can be used to evaluate the moments of the displacement. For the average we have

[TABLE]

where $\langle\Delta x|f\rangle$ denotes the average displacement for a given $f$ as in eq. (9). Making use of the expression in (16) we obtain:

[TABLE]

Similarly, for the variance one can write:

[TABLE]

which embodies the same result as in (1). In moving from the second line to the third one we have used the expression for the variance of the time spent in $+$ given in (17). As anticipated, the total dispersion involves thermal fluctuations $T\bar{\mu}$ and a term related to the mobility switching which is proportional to the square of the applied force and the variance of the time spent in a state. For the third cumulant one similarly finds

[TABLE]

The procedure extends in a simple way to higher moments.

The case of several internal states can be treated by the presented scaled cumulant generating function but requires a slight reformulation that we detail in Appendix A. The main difference is that one does not directly compute the scaled cumulant generating function $\lambda_{0}(s)$ but evaluates the various moments by means of perturbation theory expanding for small $s$.

V Multiple-scales technique

An alternative way to compute the first two moments of the displacement is to use a multiple-scale technique as done, for instance, for rotating Brownian particles in aurell2016 ; marino2016b . To this aim it is convenient to rewrite the Fokker-Planck equation (II) as

[TABLE]

where

[TABLE]

and $\mathsf{K}$ is defined in Eq. (3). In the following we shall make use of the eigendecomposition of $\mathsf{K}$ which we detail in Appendix B. Since we are seeking a diffusive effective dynamics we will introduce two time scales and adopt the the following scaling in space

[TABLE]

and require that

[TABLE]

The first time $\theta$ in (35b) is the scale of the jump process, on the order of $k^{-1}$ (jump rates in either direction, provided they are not widely different). In that time the particle on the average moves a distance about $k^{-1}\mu F$ with scatter on the order $\sqrt{T\mu k^{-1}}$. These lengths are the first spatial scale in (35a). The second spatial scale in (35a) is where we observe the motion, and the third time $\tau$ in (35b) is the time scale on which we expect to see diffusive behavior on that scale. The second time $\vartheta$ in (35b) is the time it takes for the particle to traverse the large spatial scale moving steadily with the mean velocity $\mu F$; in the present analysis this will be only needed in an intermediate technical step, and does not by itself generate new physical effects. With these definitions, $X$ is the the large scale space variable which is of order one only for very large $x$. We are interested in finding the effective dynamics in $X$. The small scale variable $\tilde{x}$ essentially corresponds to the original variables $x$ but is restricted to small scales by imposing periodic boundary conditions for $P$ over the typical length covered during the relaxation time of $a$. As a consequence of the scaling (35) and (36), the time and spatial derivatives in turn into

[TABLE]

while the transition matrix $\mathsf{K}$ remains unchanged. We treat $\varepsilon$ as a small perturbative parameter and expand $P$ in powers of $\varepsilon$,

[TABLE]

where all $P^{(i)}$ depends on the various variables as in eq. (36) . In these variables, $P^{(0)}$ is normalized to one, while all other $P^{(i)}$ with $i>0$ are normalized to zero. Plugging (37) and (38) into (30), and collecting terms of equal powers in $\varepsilon$ in the resulting expression, we obtain a hierarchy of inhomogeneous Fokker-Planck like equations of which we list the first three (order $\varepsilon^{0}$, $\varepsilon^{1}$ and $\varepsilon^{2}$):

[TABLE]

In (39), we introduced the tilde over the operators $\tilde{{\cal{L}}}^{\dagger}$ to indicate that it acts on the small scale variables $\tilde{x}$. Note that (39a) is the same equation as (30), with the essential difference, however, that $P^{(0)}$ obeys periodic boundary conditions in the variables $\tilde{x}$. We follow the standard procedure detailed, e.g., in pavliotis08 . We are interested in solutions of (39) which are stationary on small scales after short-term transients have died out. Hence, the desired solutions do not depend on $\theta$ such that we can set $\partial P^{(i)}/\partial\theta=0$ for all $i$.

Order $\varepsilon^{0}$.

After relaxation of the fastest timescale, the periodic solution in $\tilde{x}$ is given by a constant in $\tilde{x}$: $\rho(\vartheta,\tau,X)$. For the $a$ variable the system relaxes to the steady state (4) so that we can write:

[TABLE]

Order $\varepsilon^{1}$.

After relaxation of the fastest time variable the LHS of equation 39f reads: $-\left(\tilde{{\cal{L}}}^{\dagger}+\mathsf{K}\right)P^{(1)}$. For a solution to exist we need to require the RHS to be orthogonal to the left null space of the operator. Such null space is spanned by constants multiplied by the left null space of $\mathsf{K}$ which is given by $\hat{\boldsymbol{w}}=(1,\;1)$ (see Appendix B). Multiplying the RHS of equation 39f corresponds to summing its rows and the solvability condition gives:

[TABLE]

where

[TABLE]

is the same as in (IV) since $\overline{f}=w_{+}$, and where we have used that $P^{(0)}$ is independent of $\tilde{x}$. We can now plug this expression in eq. (39f) which, after the initial relaxation now reads:

[TABLE]

where, in moving to the last line we have exploited that the parameters related to the jump process ($k_{-}$, $k_{+}$, $w_{-}$, $w_{+}$, $\mu_{-}$, $\mu_{+}$, …) are independent of space and where

[TABLE]

is the right eigenvector associated with the non zero eigenvalue of $\mathsf{K}$, which as required by the solvability condition, is orthogonal to $\hat{\boldsymbol{w}}$ (see Appendix B). Since the RHS does not depend on the small spatial scale $\tilde{x}$ we make the ansatz that also $P^{(1)}$ is independent of them. Then, $P^{(1)}$ is found by applying $\mathsf{G}$ the Green’s function of the operator $\mathsf{K}$ given in eq. (89) to the RHS. The Green’s function is defined so that $\mathsf{K}\mathsf{G}=-\delta_{i}^{j}+\boldsymbol{w}\hat{\boldsymbol{w}}$. Making again use of the fact that the switching rates (and consequently $\mathsf{G}$ and $\boldsymbol{w}$) are independent of space, we can write

[TABLE]

Order $\varepsilon^{2}$.

For equation (39i) to admit a solution, we need to require its RHS to be orthogonal to the left null space of $\mathsf{K}$ spanned by $\hat{\boldsymbol{w}}$. Plugging the solution for $P^{(1)}$ (eq. (V))) and the condition (41) into the solvability condition leads to:

[TABLE]

where we have made use of the fact that $\left(\mu_{+}-\overline{\mu},\,\mu_{-}-\overline{\mu}\right)=-(\mu_{+}-\mu_{-})\frac{k_{-}}{k_{+}+k_{-}}\hat{\boldsymbol{n}}_{1}$. The connection with the theory of Taylor dispersion can be made explicit by considering the first term on the RHS of the first row of (V) and the explicit expression of $P^{(1)}$ in (V). The Taylor contribution reads:

[TABLE]

which coincides with the known expression (see e. g. Eq. (19) in Ref. kahlen2017 ) if one identifies $\mu F$ with the state dependent velocity $u$ in kahlen2017 . Composing the two solvability conditions and returning to the original variables we have the following SDE for the effective dynamics

[TABLE]

showing that the effective equation evolves with average

[TABLE]

and variance

[TABLE]

in accordance with Eqs. (IV) and (IV) obtained by using the long time limit of the cumulants of the residence time in state $+$.

VI Generalizations: space-dependent switching rates and several internal states

In many systems of interest, (e.g., molecular motors watanabe2012mechanical ; pietzonka2014fine ; kawaguchi2014nonequilibrium ) the transition rates between the internal states may depend on the particle position. In general, this may greatly complicate the treatment. However, under certain assumptions on the spatial dependence, the multiple scale method presented in previous section can be applied. For instance, if the switching rates vary appreciably only on the large spatial scales on which the effective diffusion is observed, the operator $\mathsf{K}$ will depend parametrically on the large spatial coordinate $X$ but will act on the internal state in the same way. The zeroth-order equation (40) then still applies, the only difference being that $\boldsymbol{w}$ is a function of $X$. The first order solvability equation (41) is also unchanged which means that the drift velocity is still $F\bar{\mu}$, where $\bar{\mu}$ (which however now depends on $X$) is given in (42). The equation to solve on first order (45) is also formally unchanged, but its coefficients now depend on $X$ and so does the Green’s function $G$.

The first point where this type of large spatial scale dependence of the switching rates enters the analysis is therefore (V) which instead will read

[TABLE]

The Green’s function can be brought inside the derivative in (60), at the price of a correlation term proportional $\frac{\partial\mathsf{G}}{\partial X}$. The upshot is hence that on the large scale motion is diffusive with formally the same diffusion coefficient as in (55) (now space dependent and following the Itô prescription), but with a correction to the drift velocity reading:

[TABLE]

In principle, a solution with a multiple-scale approach can be found also if the rates depend periodically on space with a period much smaller than the distance travelled by the particle between different jumps in the mobility state. The case in which the rates vary on the same length scale at which the particle moves is more challenging and may not be solvable analytically with a multiple scale approach. An extended discussion of a related problem concerning the case of rotational and translational diffusion in periodic potentials can be found in marino2016dynamics (see paper IV therein).

The multiple-scale analysis can be straightforwardly generalized to the case in which there are more than two discrete states. In this Section we extend our formalism to the case of $N$ internal states. The effective displacement reads as Eq. (56) with the average mobility now displaying the sum over all states: $\overline{\mu}=\sum_{i}\mu_{i}w_{i}$ where $w_{i}$ represents the steady state probability of the discrete process being in state $i$. Similarly, the first contribution to the variance, remains proportional to the average mobility as in Eq. (57). The Taylor dispersion contribution to the variance maintains the structure of Eq. (54) and can be written as

[TABLE]

where we have introduced the matrix $\tilde{\mathsf{\mu}}=\delta_{i}^{j}(\mu_{i}-\overline{\mu})$ and, as before, $\mathsf{G}$ is the Green’s function of the transition matrix $\mathsf{K}$ (on the subspace on non-zero eigenvalues) and $\boldsymbol{w}$ and $\hat{\boldsymbol{w}}$ are respectively the steady state solution and a unit vector but now extended to $N$ states. The main difference is a computational one since explicitly finding the Green’s function can become complicated for larger matrices. An analysis by scaled cumulant function for several (discrete or continuous) internal states is carried out in Appendix A and Appendix C.

VII Coarse-grained properties as a measurement device

As mentioned in the introduction, the fact that the particle dynamics is influenced by the properties of the friction switching makes it possible to exploit measurements of the particle position to gain insight about the internal dynamics and evaluate its transition rates. Remarkably, despite the fact that we are observing the particle on time scales that are longer than the relaxation time of the internal process, we do not only have access to the equilibrium properties of the internal variable but also to its dynamical ones. This is possible because we are driving the system away from equilibrium by applying the constant force $F$. More precisely, from the long-time average displacement of the particle (IV) it is possible to measure the average time spent in a state, i. e., the equilibrium probability of the internal process $w_{+}$ and $w_{-}$ which involve the ratio of the transition rates (4). To perform such measurement it is necessary to know the friction values corresponding to the two states $\mu_{+}$ and $\mu_{-}$ and the externally applied force $F$. The long-term particle dispersion can be exploited to evaluate the magnitude of the friction transition rates. Expression (1) is the sum of two contributions: a first one being the average of the diffusion constant in the two states (depending on the rate ratio) and another one caused by the combination of the external force and the friction fluctuations which depends on the variance of the time spent in a given state and involves the magnitude of the transition rates. The relative importance of the switching dependent contribution depends on how different the two friction coefficients are and on how similar the time scales of diffusion and reaction are. Let us introduce the time on which diffusion and drift have comparable effects as $\tau_{D}=\overline{\mu}T/(\overline{\mu}F)^{2}$ and the typical time of the friction switching: $\tau_{R}=\lambda_{1}^{-1}=(k_{+}+k_{-})^{-1}$. We then see that the ratio between the fluctuations due to mobility switching and the thermal ones is proportional to the ratio of the two time scales:

[TABLE]

Such dependence on the time scale can be traced back to the fact that if the friction dynamics is faster than the diffusive one, many transitions will occur before any considerable particle displacement takes place. This averages out fluctuations in the time spent in a given friction state and consequently reduces variance. In such case, the particle dynamics will follow Eq. (2) with the state dependent mobility $\mu_{a}$ replaced by the average one $\overline{\mu}$ in analogy with the cases discussed in pavliotis08 ; bo2017 (see pietzonka2016extreme for a detailed discussion). Let us consider a concrete example with a particle, immersed in water at room temperature, switching between two spherical conformations respectively with radius $0.6\,\mathrm{\mu m}$ and $0.4\,\mathrm{\mu m}$ with identical rates $k$, subject to a force of $F=100\,\mathrm{fN}$. On average it will move by about $11\,\mathrm{\mu m\;s^{-1}}$. If the switching occurs at a rate $k=10\,\mathrm{s^{-1}}$ the position variance will grow as $1.4\,\mathrm{\mu m^{2}\;s^{-1}}$ with the force dependent contribution making up about $36\%$ of the total.

As we have seen one can derive the moments of the particle displacements also for cases with more than two internal states. In this case, however, knowing the average and variance of the displacement will no longer be sufficient to estimate the transition rates since there are more unknown rates than known moments. One then has to proceed to measure (and compute) higher moments, a procedure which will cease to be practical as the number of states increases. However, for limited systems this may be possible. Consider for instance, a system composed of three states $+\rightleftharpoons 0\rightleftharpoons-$ with four independent rates. With the presented method, their estimation then requires four moments. Reliably measuring higher moments can be experimentally challenging but moments up to the fourth moment can be obtained. For instance, the ratio of the fourth to second moment is used to estimate the non Gaussianity of anomalous diffusion (see e.g. Ref. metzler2016non and references therein).

VIII Conclusions and discussion

In this work we have reconsidered a result dating back more than half a century, that the diffusion of a particle with mobility dependent on an internal state has two conceptually quite different components. The first term in the effective diffusion is of the standard thermal type and satisfies an Einstein relation with the effective mobility. The second term is on the other hand quadratic in the externally applied force and depends differently on the kinetic rates of transitions between the internal states. We have pointed out that measuring the first two moments of displacement hence determines both the rates in both directions when there are only two internal states and discussed its possible use as an in vitro measurement device from a modern perspective. We have generalized these results to more than two internal states, and, in Appendices, to a continuum of internal states.

The results of this paper have been derived by a scaled cumulant expansion (large deviations) and by a multiple-scales calculation. The extension to the continuous case enabled us to show that the first and second order moments of translation of a rotating Brownian particle can be derived by the scaled cumulant expansion. The continuous internal state is then the orientation of the particle, an element of the three-dimensional rotation group SO(3). In a previous contribution aurell2016 we derived this result by a multiple-scales approach. A conclusion of this work is thus that for the problems considered the multiple-scales approach can be embedded in the wider mathematical framework of large deviation theory.

Acknowledgements.

We thank Raffaele Marino, Ralf Eichhorn and Erwin Frey for stimulating discussions and Udo Seifert and Ramin Golestanian for constructive remarks. This research was supported by the Academy of Finland through its Center of Excellence COIN (grant no. 251170). SB wishes to thank Peking University for hospitality.

Appendix A Generalization of the scaled cumulant generating function approach to several internal states

The presence of several internal states requires a slight extension of the scaled cumulant generating function approached presented in Section IV. The starting point is unchanged, as we have that the particle evolution at a fixed value of the internal variable is a Gaussian as given in Eq. (6). The probability of observing a displacement conditioned on a trajectory of the internal variable now depends on the fraction of time spent in each state. It is convenient to introduce the column vector $\boldsymbol{f}$ with entries $f_{i}$ given by the fraction of time spent in each state $i$ and obviously normalized $\sum_{i=1}^{N}f_{i}=1$. The expression for the probability is the same of the one in Eq. (8) for the two states system as well as Eqs. (9):

[TABLE]

and (10) with the notable difference that $f$ is replaced by $\boldsymbol{f}$ and that

[TABLE]

The main difference with the simpler case is that one needs to perform averages over the joint distribution $P(\boldsymbol{f})=P(f_{1},\,f_{2},\;\ldots)$. The average displacement, in analogy with Eq. (IV), reads:

[TABLE]

where the overline refers to the average over the joint probability $P(\boldsymbol{f})$. In the long-time limit we have that the average fraction of time spent in a state coincided with the probability of being in that state $\overline{f_{i}}=w_{i}$. For the variance of the displacement one has

[TABLE]

requiring the evaluation of the covariance matrix of the fraction of time spent in a state

[TABLE]

To this aim, we can make use of the generating function (in analogy to eq.(13))

[TABLE]

and exploit the fact that

[TABLE]

For long times this can be solved by finding the leading eigenvalue $\lambda_{0}(\boldsymbol{s})$ of the matrix:

[TABLE]

where we have introduced

[TABLE]

This can be achieved by means of perturbation theory

[TABLE]

where $\lambda_{k}(\boldsymbol{0})$ are the eigenvalues of the unperturbed matrix $\mathsf{K}$ and $\boldsymbol{m}_{k}$, $\hat{\boldsymbol{n}}_{k}$ respectively its right and left eigenvectors. For consistence with earlier notation we have denoted the eigenvectors of the zero eigenvector as $\boldsymbol{m}_{0}=\boldsymbol{w}$, $\hat{\boldsymbol{n}}_{0}=\hat{\boldsymbol{w}}$. We have also introduced the Green’s function $-\sum_{k>0}\frac{1}{\lambda_{k}(\boldsymbol{0})}\boldsymbol{m}_{k}\hat{\boldsymbol{n}}_{k}=\mathsf{G}$. The covariance then reads

[TABLE]

where no summation is to be taken over the repeated indices. With this expression the term proportional to the force squared in (A), reads

[TABLE]

corresponding with the diffusion related to the term in Eq. (65), where we have exploited the biorthogonality $\overline{\mu}^{2}\hat{\boldsymbol{w}}\mathsf{G}\boldsymbol{w}=0$.

Appendix B Eigen decomposition of $\mathsf{K}$ for two states

In this section we will present the eigen-decomposition of the matrix $\mathsf{K}$. Since it generates a Markov process, the largest eigenvalue of $\mathsf{K}$ is $\lambda_{0}=0$ and its associate right eigenvector gives the steady state distribution

[TABLE]

and its left eigenvector

[TABLE]

ensures probability conservation. With

[TABLE]

and

[TABLE]

The second eigenvalue is given by

[TABLE]

and its eigenvectors are

[TABLE]

where

[TABLE]

and

[TABLE]

This provides a biorthonormal decomposition with:

[TABLE]

and the completeness relation

[TABLE]

The Green’s function

[TABLE]

ensures that for vectors orthogonal to the left null space of $\mathsf{K}$ (parallel to $\boldsymbol{m}_{1}$)

[TABLE]

Appendix C Large deviation analysis with internal continuous degree of freedom

In this Appendix we make the connection to our earlier work on coarse-graining of diffusion in space and over orientations aurell2016 ; marino2016b . The starting point is the same as in the preceding discussion, except that the internal variable $\alpha$ representing the internal state now lives in a compact manifold ${\cal M}$ with volume element $\sqrt{g}$. We also assume that the spatial motion of the particle takes place in $D$-dimensional space. Rotations in $D=2$ can hence be represented by an internal coordinate on a circle (Lie group $SO(2)$), rotations in $D=3$ by an internal coordinate in the solid upper hemisphere (Lie group $SO(3)$), and so on. In marino2016b the case of motion in two-dimensional complex space ($D=4$) and the group $SU(2)$ was considered as well. Whatever the internal state the process is described by

[TABLE]

As in the previous discussion $\mathbf{\mu}(\alpha)$ in (91) is a mobility matrix and $\mathbf{\mu}^{\frac{1}{2}}$ is its matrix square root, $\text{d}W_{j}$ is the standard Wiener increment, $v_{i}(\alpha)=\mathbf{\mu}(\alpha)_{ij}F_{i}$ is a drift velocity and $F$ is an applied force, assumed constant in space and time. We assume that $\mathbf{\mu}(\alpha)$ depends smoothly on $\alpha$.

Now suppose that $\alpha$ obeys a dynamic law independent of spatial position $x$ such that in a total time $t$ the coordinate $\alpha$ can be found in a set $A\subset{\cal M}$ for a time $t_{A}$. In parallel with the discussion in the main text we can introduce $f_{A}(t)=\frac{t_{A}}{t}$, the fraction of the time $\alpha$ can be found in $A$ in a total time $t$. Suppose further that ${\cal M}$ is divided up in a collection of small sets $A_{k}$ with zero-measure intersections and which together cover ${\cal M}$. Approximately the motion must then be same as with a finite number of internal states discussed in the preceding Section A.

In the limit that all the sets $A_{k}$ are small we assume that there is a normalized empiric measure $m(\alpha)\sqrt{g}$ smoothly related to the volume element such that $f_{A_{k}}=\int_{A_{k}}\sqrt{g}m(\alpha)$ and $\int_{\cal M}\sqrt{g}m(\alpha)=1$. In analogy with the previous case this leads to an overall mobility dependent on the empiric measure

[TABLE]

and an overall drift velocity

[TABLE]

In both (92) and (93) the integrals are performed with the empirical measure, and $\mathbf{\mu}[m]$ and $\boldsymbol{v}[m]_{i}$ are hence both linear functionals of the scalar function $m(\alpha)$.

We now further postulate a functional $Q[m(\alpha),t]$ such that ${\cal D}mQ[m]$ is the probability of an empirical measure in a functional volume element ${\cal D}m$ around $m(\alpha)\sqrt{g}$. $Q$ could be taken the limit of finite-dimensional probability distributions $Q[\{f_{V_{k}}\},t]$ where $\{f_{V_{k}}\}$ is a finite set of empirical frequencies analogous to $\boldsymbol{f}$ in the previous section. We define the average of the empirical measures as

[TABLE]

The effective velocity will then be the average with respect to $Q$ of the drift velocity in (93) i.e.

[TABLE]

In the long time limit $\overline{m(\alpha)}$ will be the same as the steady state distribution $\rho^{ss}(\alpha)$, as will be verified below.

Introducing for compactness $\delta\boldsymbol{X}=\Delta\boldsymbol{X}-t\boldsymbol{v}^{eff}$ we can compute the centered second moment

[TABLE]

The two functional integrals in (C) can respectively be written

[TABLE]

and

[TABLE]

We introduce the moment generating functional of the time spent at internal position $\alpha$, conditional on ending up at internal position $\alpha^{\prime}$ at time $t$. This moment generating functional depends on a function $s(\alpha)$ and is defined as

[TABLE]

We introduce also the probability $Q[m]^{\alpha^{\prime}}$ of observing a frequency distribution in a functional volume element around measure $m(\alpha)$ for a process ending at $\alpha^{\prime}$ at time $t$. Then

[TABLE]

We now assume that the time development of $\alpha$ is a stochastic process with generator ${\cal L}$. The moment generating functional then obeys

[TABLE]

where ${\cal L}^{\dagger}$ is the adjoint of the generator (Fokker-Planck operator for a diffusion process) acting at internal coordinate $\alpha^{\prime}$. For long times the moment generating function will be almost independent of $\alpha^{\prime}$ and dominated by $t\lambda_{0}[s]$, where $\lambda_{0}[s]$ is the largest eigenvalue of the operator ${\cal L}^{\dagger}+s(\alpha)$. If $s=0$ this largest eigenvalue must be zero, corresponding to the stationary probability distribution for the stochastic process in $\alpha$. Since ${\cal M}$ is assumed compact this eigenvalue must further be separated by a gap from all the decaying states, and $\lambda_{0}[s]$ must therefore be close to zero for $s$ sufficiently small.

To evaluate the terms in (98) and (99) we consider as before in the long term limit of the generating functional and its functional logarithmic derivatives evaluated at $s=0$:

[TABLE]

These functional derivatives are in the long time limit as $t\frac{\delta\lambda[s]}{\delta s(\alpha)}$ and $t\frac{\delta^{2}\lambda[s]}{\delta s(\alpha)\delta s(\alpha^{\prime})}$ and the eigenvalue can be computed in in time-independent perturbation theory

[TABLE]

where $m_{k}$ and $n_{k}$ are the right and left eigenfunctions of ${\cal L}^{\dagger}$ corresponding to eigenvalue $\lambda_{k}[0]$ and the scalar product is $(f,g)=\int_{\cal M}\sqrt{g}\,fg$. The average of the empirical measure is therefore (in the long-time limit) the first functional derivative of (C) with respect to $s(\alpha)$, which gives

[TABLE]

The left eigenfunction $n_{0}$ is the zero mode of the generator ${\cal L}$. For diffusions, where all the derivatives in ${\cal L}$ are on the right, $n_{0}$ must then be a constant. By orthonormality $(n_{0},m_{0})=1$ and it is convenient to take $n_{0}=1$ and $m_{0}$ a normalized probability distribution on ${\cal M}$ i.e. $\int_{\cal M}\sqrt{g}\,m_{0}=1$. This normalized probability distribution is at the same time the right eigenvector of the adjoint of the generator with eigenvalue zero i.e. the stationary probability distribution $\rho^{ss}(\alpha)$ of the stochastic process, as already observed above. For the first term in (C) we combine (98) and (104) and find a first term in an effective diffusion tensor

[TABLE]

where the effective mobility tensor has been introduced in (98)

The covariance of the empirical measure is

[TABLE]

where we have introduced a Green’s function

[TABLE]

and used that $n_{0}(\alpha)$ is a constant that we can set to one. The Green’s function satisfies the equation

[TABLE]

where $\Pi_{0}$ is a projection operator that acts on a function $f$ as $[\Pi_{0}f](\alpha)=(n_{0},f)m_{0}(\alpha)$.

The second term in (C) i.e. the contribution from the functional integral in (99) involves two mobility matrices $\mu_{ik}(\alpha)$ and $\mu_{jl}(\alpha^{\prime})$ and the Green’s function ${\cal G}(\alpha,\alpha^{\prime})$. A convenient representation is to first introduce an auxiliary function

[TABLE]

Acting with ${\cal L}^{\dagger}$ on $\lambda$ gives

[TABLE]

We observe that

[TABLE]

where the inner integral is

[TABLE]

and the last equality follows because all the right eigenfunctions $m_{k}$ with $k>0$ must be orthogonal to $n_{0}$ which is a constant. The second term can thus be written

[TABLE]

where the fourth-order tensor $A$ is

[TABLE]

For the cases studied in aurell2016 ; marino2016b $\mathbf{\mu}_{jl}(\alpha)$ equals $[R(\alpha)\mu^{B}R^{-1}(\alpha)]_{jl}$ where $R(\alpha)$ is a rotation matrix and $\mu^{B}$ is a mobility matrix in the body frame, and $m_{0}$ is the uniform (Haar) measure. As discussed therein $\mathbf{\mu}^{eff}_{jl}$ then equals $\frac{1}{3}\hbox{Tr}[\mathbf{\mu}]\mathbf{1}_{jl}$ and the quantity in the parenthesis on the right-hand side of (110) is the traceless part of the mobility matrix. In aurell2016 and marino2016b this quantity was denoted $\tilde{\gamma}^{-1}$. Eqs. (110), (112) and (113) generalize the corresponding expressions from aurell2016 to cases where steady state probability is not necessarily uniform.

Appendix D Higher moments for continuous internal variables

In this Appendix we compute the centered third moment with the scaled cumulant expansion technique and show that it is linear in time $t$. It can be written as a sum of three terms given below in (116) and (118). The starting point is

[TABLE]

where the second term generalizes (106) and (99) and can be written

[TABLE]

Comparing to (99) (112) and (113) one sees that the first term in (114) is linear in time and given by

[TABLE]

The tensor $A$ is given in (113) above. Turning to the other term, the expression inside the inner square bracket of (118) is $\frac{\delta^{3}\log G[s]}{\delta s(\alpha)\delta s(\alpha^{\prime})\delta s(\alpha^{\prime\prime})}$ which is equal to $t\frac{\delta^{3}\lambda_{0}[s]}{\delta s(\alpha)\delta s(\alpha^{\prime})\delta s(\alpha^{\prime\prime})}$. Also this term is therefore linear in time. To compute it we thus need the third order perturbation of an energy which is

[TABLE]

where ${\cal G}^{(2)}(\alpha^{\prime},\alpha^{\prime\prime})=\sum_{n>0}m_{n}(\alpha^{\prime})n_{n}(\alpha^{\prime\prime})\frac{1}{(\lambda_{n}[0])^{2}}$ satisfies $\left({\cal L}^{\dagger}\right)^{2}{\cal G}^{(2)}(\alpha^{\prime},\alpha^{\prime\prime})=\delta(\alpha^{\prime}-\alpha^{\prime\prime})-\Pi_{0}$. This means that the the second term in (114) can be written as a sum of two terms

[TABLE]

where the two sixth order tensors are given by

[TABLE]

and

[TABLE]

For the term $B^{(1)}$ we can use the same auxiliary quantity as above

[TABLE]

and further another auxiliary quantity

[TABLE]

which satisfies $({\cal L}_{0}\lambda^{L}_{il})(\alpha^{\prime})=\mu_{il}(\alpha^{\prime})-\mu_{il}^{eff}$ and which together give

[TABLE]

For the term $B^{(2)}$ we introduce yet another auxiliary quantity

[TABLE]

in terms of which we have

[TABLE]

When $m_{0}$ is the Haar measure the averages in (123) (125) can be computed by invariance arguments as in aurell2016 and marino2016b .