100 years after Smoluchowski: stochastic processes in cell biology

Obviously, the infinite velocities of the Brownian trajectories contradict physics. Specifically, according to the Waterston–Maxwell equipartition theorem [20], the root mean square (RMS) velocity $\bar{v}=\sqrt{\mathbb{E}{{v}^{2}}}$ of a suspended particle should be determined by the equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{m}{2}{{{\bar{v}}}^{2}}=\frac{3kT}{2}. \end{array}\end{eqnarray} \tag{ 7 }$$

Each component of the velocity vector has the same variance, so that

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{m}{2}\bar{v}_{x,y,z}^{2}=\frac{kT}{2}, \end{array}\end{eqnarray} \tag{ 8 }$$

which is the one-dimensional version of equation (7). The RMS velocity comes out to be about 8.6 cm s⁻¹ for the particles used in Svedberg's experiment [22]. Einstein argued in 1907 and 1908 [1] that there is no possibility of observing this velocity, because of the very rapid viscous damping, which can be calculated from the Stokes formula. The velocity of such a particle would drop to 1/10 of its initial value in about $3.3\times {{10}^{-7}}$ s. Therefore, Einstein argued, in the period τ between observations the particle must get new impulses to movement by some process that is the inverse of viscosity, so that it retains a velocity whose RMS average is $\bar{v}$. Between consecutive observations these impulses alter the magnitude and direction of the velocity in an irregular manner, even in the extraordinarily short time of $3.3\times {{10}^{-7}}$ s. According to this theory, the RMS velocity in the interval τ has to be inversely proportional to $\sqrt{\tau}$; that is, it increases without limit as the time interval between observations becomes smaller.

In 1908 Langevin [2] offered an alternative approach to the problem of the Brownian motion. He assumed that the dynamics of a free Brownian particle is governed by the frictional force $-6\pi a\eta v$ and by a fluctuational force $\Xi$ that results from the random collisions of the Brownian particle with the molecules of the surrounding fluid, after the frictional force is subtracted. This force is random and assumes positive and negative values with equal probabilities. It follows that Newton's second law of motion for the Brownian particle is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} m\ddot{{x}}\,=-6\pi a\eta \overset{\centerdot}{{x}}\,+ \Xi. \end{array}\end{eqnarray} \tag{ 9 }$$

Denoting $v=\overset{\centerdot}{{x}}\,$ and multiplying equation (9) by x, we obtain

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{m}{2}\frac{{{\text{d}}^{2}}}{\text{d}{{t}^{2}}}{{x}^{2}}-m{{v}^{2}}=-3\pi a\eta \frac{\text{d}}{\text{d}t}{{x}^{2}}+ \Xi x. \end{array}\end{eqnarray} \tag{ 10 }$$

Averaging under the assumption that the fluctuational force $\Xi$ and the displacement of the particle x are independent, we obtain

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{m}{2}\frac{{{\text{d}}^{2}}}{\text{d}{{t}^{2}}}\mathbb{E}{{x}^{2}}+3\pi a\eta \frac{\text{d}}{\text{d}t}\mathbb{E}{{x}^{2}}=kT, \end{array}\end{eqnarray} \tag{ 11 }$$

where (8) has been used. The solution is given by $\text{d}\mathbb{E}{{x}^{2}}/\text{d}t=kT/3\pi a\eta +C{{\text{e}}^{-6\pi a\eta t/m}}$, where C is a constant. The time constant in the exponent is ${{10}^{-8}}\,$ sec, so the mean square speed decays on a time scale much shorter than that of observations. It follows that $\mathbb{E}{{x}^{2}}-\mathbb{E}x_{0}^{2}=\left(kT/3\pi a\eta \right)t$. This, in turn (see (4)), implies that the diffusion coefficient is given by $D=kT/6\pi a\eta$, as in Einstein's equation (1).

Langevin's equation (9) is a stochastic differential equation, because it is driven by a random force $\Xi$. If additional fields of force act on the diffusing particles (e.g. electrostatic, magnetic, gravitational, etc), Langevin's equation is modified to include the external force, F(x,t), say, [12], [11],

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} m\ddot{{x}}\,+ \Gamma \overset{\centerdot}{{x}}\,-F(x,t)= \Xi, \end{array}\end{eqnarray} \tag{ 12 }$$

where $\Gamma =6\pi a\eta$ is the friction coefficient of a diffusing particle. We denote the dynamical friction coefficient (per unit mass) $\gamma = \Gamma /m$. If the force can be derived from a potential, $F=-\nabla U(x)$, Langevin's equation takes the form

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} m\ddot{{x}}\,+ \Gamma \overset{\centerdot}{{x}}\,+\nabla U(x)= \Xi. \end{array}\end{eqnarray} \tag{ 13 }$$

The main mathematical difference between the two approaches is that Einstein assumes that the displacements $\Delta$ are independent, whereas Langevin assumes that the random force $\Xi$ and the displacement x are independent. The two theories are reconciled in section 3.3 below.

To investigate the statistical properties of the fluctuating force $\Xi$, Langevin made the following assumptions.

• (i)  
  The fluctuating force $\Xi$ is independent of the velocity v.
• (ii)  
  $\Xi$ changes much faster than v.
• (iii)  
  $\langle \Xi \rangle =0$.
• (iv)  
  The accelerations imparted in disjoint time intervals $\Delta {{t}_{1}}$ and $\Delta {{t}_{2}}$ are independent.

These conditions define the noise $\Xi (t)$ (so called white noise) as the nonexistent derivative of Einstein's Brownian motion X(t). The conditional probability distribution function of the velocity process of a Brownian particle (PDF), given that it started with velocity v₀ at time t  =  0, is defined as $P\left(v,t\,|\,{{v}_{0}}\right)=\text{Pr}\left\{v(t)<v\,|\,{{v}_{0}}\right\}$ and the conditional probability density function is defined by

$$\begin{eqnarray}&&p\,\left(v,t\,|\,{{v}_{0}}\right)=\frac{\partial P\left(v,t\,|\,{{v}_{0}}\right)}{\partial v}.\end{eqnarray}$$

In higher dimensions, we denote the displacement vector $\boldsymbol{x}={{\left({{x}_{1}},{{x}_{2}},\ldots,{{x}_{d}}\right)}^{T}}$, the velocity vector $\overset{\centerdot}{{\boldsymbol{x}}}\,=\boldsymbol{v}={{\left({{v}_{1}},{{v}_{2}},\ldots,{{v}_{d}}\right)}^{T}}$, the random force vector $\boldsymbol{\Xi}={{\left({{ \Xi }_{1}},{{ \Xi }_{2}},\ldots,{{ \Xi }_{d}}\right)}^{T}}$, the PDF

$$\begin{eqnarray}&&P\left(\boldsymbol{v},t\,|\,\boldsymbol{v}_{{0}}\right)=\text{Pr}\left\{{{v}_{1}}(t)<{{v}_{1}},{{v}_{2}}(t)<{{v}_{2}},\ldots,{{v}_{d}}(t)<{{v}_{d}}\,|\,\boldsymbol{v}(0)=\boldsymbol{v}_{{0}}\right\},\end{eqnarray}$$

and the probability density function (pdf)

$$\begin{eqnarray}&&p\,\left(\boldsymbol{v},t\,|\,\boldsymbol{v}_{{0}}\right)=\frac{{{\partial}^{n}}P\left(\boldsymbol{v},t\,|\,\boldsymbol{v}_{{0}}\right)}{\partial {{v}_{1}}\partial {{v}_{2}},\ldots \partial {{v}_{d}}}.\end{eqnarray}$$

The conditioning implies that the initial condition for the pdf is $p\,\left(\boldsymbol{v},t\,|\,\boldsymbol{v}_{{0}}\right)\to \delta \left(\boldsymbol{v}-\boldsymbol{v}_{{0}}\right)$ as $t\to 0$. According to the Waterston–Maxwell theory, when the system is in thermal equilibrium, the velocities of free Brownian particles have the Maxwell–Boltzmann pdf; that is,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{t\to \infty}{{\lim}}\,p\,\left(\boldsymbol{v},t\,|\,\boldsymbol{v}_{{0}}\right)={{\left(\frac{m}{2\pi kT}\right)}^{3/2}}\exp \left\{-\frac{m|\boldsymbol{v}{{|}^{2}}}{2kT}\right\}. \end{array}\end{eqnarray} \tag{ 14 }$$

The solution of the Langevin equation (9) for a free Brownian particle is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{v}(t)=\boldsymbol{v}_{{0}}{{\text{e}}^{-\gamma t}}+\frac{1}{m}{\int}_{0}^{t}{{\text{e}}^{-\gamma (t-s)}}\boldsymbol{\Xi}(s)\,\text{d}s. \end{array}\end{eqnarray} \tag{ 15 }$$

The integral (15) makes sense, if it is integrated once by parts. However, if the factor multiplying $\boldsymbol{\Xi}$ in (15) is non-differentiable as well, e.g. if it is $\boldsymbol{x}(t)$, integration by parts is insufficient to make sense of the stochastic integral and a more sophisticated definition is needed. Such a definition was given by Itô in 1944 [29] (see also [18]).

To interpret the stochastic integral in the one-dimensional (15), we make a short mathematical digression on the definition of integrals of the type ${\int}_{0}^{t}g(s) \Xi (s)\,\text{d}s$, where g(s) is a deterministic integrable function. Such an integral is defined as the limit of finite Riemann sums of the form

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {\int}_{0}^{t}g(s) \Xi (s)\,\text{d}s=\underset{ \Delta {{s}_{i}}\to 0}{{\lim}}\,\underset{i}{\sum}\,g\left({{s}_{i}}\right) \Xi \left({{s}_{i}}\right)\, \Delta {{s}_{i}}, \end{array}\end{eqnarray} \tag{ 16 }$$

where $0={{s}_{0}}<{{s}_{1}}<\cdots <{{s}_{N}}=t$ is a partition of the interval [0,t]. According to the assumptions about $\Xi$, if we choose $\Delta {{s}_{i}}= \Delta t=t/N$ for all i, the Gaussian increments $\Delta {{b}_{i}}= \Xi \left({{s}_{i}}\right)\, \Delta {{s}_{i}}$ are independent identically distributed (i.i.d.) random variables. Einstein's observation (see the beginning of section 3) that the RMS velocity on time intervals of length $\Delta t$ are inversely proportional to $\sqrt{ \Delta t}$, implies that if the increments ${{b}_{i}}= \Xi \left({{s}_{i}}\right)\, \Delta {{s}_{i}}$ are chosen to be normally distributed, their mean must be zero and their covariance matrix must be $\langle \Delta {{b}_{i}} \Delta {{b}_{j}}\rangle =q \Delta t{{\delta}_{ij}}$ with q a parameter to be determined. We write $\Delta {{b}_{i}}\sim \mathcal{N}\left(0,q \Delta t\right)$. Then $g\left({{s}_{i}}\right) \Xi \left({{s}_{i}}\right)\, \Delta {{s}_{i}}\sim \mathcal{N}\left(0,|g\left({{s}_{i}}\right){{|}^{2}}q \Delta t\right)$, so that ${\sum}_{i}g\left({{s}_{i}}\right) \Xi \left({{s}_{i}}\right)\, \Delta {{s}_{i}}\sim \mathcal{N}\left(0,\sigma _{N}^{2}\right)$, where $\sigma _{N}^{2}={\sum}_{i}|g\left({{s}_{i}}\right){{|}^{2}}q \Delta {{s}_{i}}$. As $\Delta t\to 0$, we obtain ${{\lim}_{ \Delta t\to 0}}\sigma _{N}^{2}=q{\int}_{0}^{t}{{g}^{2}}(s)\,\text{d}s$ and ${\int}_{0}^{t}g(s) \Xi (s)\,$ $\,\text{d}s\sim \mathcal{N}\left(0,{{\sigma}^{2}}\right)$, where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\sigma}^{2}}=q{\int}_{0}^{t}{{g}^{2}}(s)\,\text{d}s. \end{array}\end{eqnarray} \tag{ 17 }$$

To interpret (15), we use (17) with $g(s)={{\text{e}}^{-\gamma (t-s)}}$ and obtain

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\sigma}^{2}}=\frac{q}{2\gamma}\left(1-{{\text{e}}^{-2\gamma t}}\right). \end{array}\end{eqnarray} \tag{ 18 }$$

Returning to the velocity vector $\boldsymbol{v}(t)$, we obtain from the above considerations

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{v}(t)-\boldsymbol{v}_{{0}}{{\text{e}}^{-\gamma t}}\sim \mathcal{N}\left(\mathbf{0},{{\sigma}^{2}}\boldsymbol{I}\right) \end{array}\end{eqnarray} \tag{ 19 }$$

with ${{\sigma}^{2}}$ given by (18). Finally, the condition (14) implies that $q=2\gamma kT/m$, so that in 3D the mean energy is as given in equation (7). Itô's construction of the stochastic integral allows g(t) to be stochastic, but independent of the Brownian increments $\Xi \left({{s}_{i}}\right) \Delta {{s}_{i}}=X\left({{s}_{i}}+ \Delta {{s}_{i}}\right)-X\left({{s}_{i}}\right)$.

In the limit $\gamma \to \infty$ the acceleration $\gamma \boldsymbol{v}(t)$ inherits the properties of the random acceleration $\boldsymbol{\Xi}(t)$ in the sense that for different times ${{t}_{2}}>{{t}_{1}}>0$ the accelerations $\gamma \boldsymbol{v}\left({{t}_{1}}\right)$ and $\gamma \boldsymbol{v}\left({{t}_{2}}\right)$ become independent. In fact, from equations (15) and $\Delta {{b}_{i}}\sim \mathcal{N}\left(0,q \Delta t\right)$, we find that

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\gamma \to \infty}{{\lim}}\,\mathbb{E}\left[\gamma \boldsymbol{v}\left({{t}_{1}}\right)\centerdot \gamma \boldsymbol{v}\left({{t}_{2}}\right)\right]=0 \end{array}\end{eqnarray}$$

and a similar result for $0<{{t}_{2}}<{{t}_{1}}$. It follows that for $\gamma \left({{t}_{1}}\wedge {{t}_{2}}\right)\gg 1$,

$$\begin{eqnarray}&&\mathbb{E}\left[\gamma \boldsymbol{v}\left({{t}_{1}}\right)\centerdot \gamma \boldsymbol{v}\left({{t}_{2}}\right)\right]=\frac{\gamma q}{{{m}^{2}}}{{\text{e}}^{-\gamma |{{t}_{2}}-{{t}_{1}}|}}\left(1+o(1)\right)=\frac{2q}{{{m}^{2}}}\delta \left({{t}_{2}}-{{t}_{1}}\right)\left(1+o(1)\right),\end{eqnarray} \tag{ 20 }$$

because for t₁  >  0

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\gamma \to \infty}{{\lim}}\,{\int}_{0}^{\infty}f\left({{t}_{2}}\right)\frac{\gamma}{{{m}^{2}}}{{\text{e}}^{-\gamma |{{t}_{2}}-{{t}_{1}}|}}\left(1+o(1)\right)\,\text{d}{{t}_{2}}=\frac{2}{{{m}^{2}}}f\left({{t}_{1}}\right) \end{array}\end{eqnarray} \tag{ 21 }$$

for all test functions f(t) in ${{\mathbb{R}}^{+}}$.

The displacement of a free Brownian particle is obtained from integration of the velocity process,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{x}(t)=\boldsymbol{x}_{{0}}+{\int}_{0}^{t}\boldsymbol{v}(s)\,\text{d}s. \end{array}\end{eqnarray} \tag{ 22 }$$

Using the expression (15) in equation (22) and changing the order of integration in the resulting iterated integral, we obtain

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{x}(t)-\boldsymbol{x}_{{0}}-\boldsymbol{v}_{{0}}\frac{1-{{\text{e}}^{-\gamma t}}}{\gamma}={\int}_{0}^{t}g(s)\boldsymbol{\Xi}(s)\,\text{d}s, \end{array}\end{eqnarray} \tag{ 23 }$$

where $g(s)=\left(1-{{\text{e}}^{-\gamma (t-s)}}\right)/m\gamma$.

Reasoning as above, we find that the stochastic integral in equation (23) is a normal variable with zero mean and covariance matrix $\boldsymbol{\Sigma}={{\sigma}^{2}}\boldsymbol{I}$, where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\sigma}^{2}}=q{\int}_{0}^{t}{{g}^{2}}(s)\,\text{d}s=\frac{q}{2{{\gamma}^{3}}}\left(2\gamma t-3+4{{\text{e}}^{-\gamma t}}-{{\text{e}}^{-2\gamma t}}\right). \end{array}\end{eqnarray} \tag{ 24 }$$

The moments of the displacement are

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathbb{E}\left[\boldsymbol{x}(t)-\boldsymbol{x}_{{0}}-\boldsymbol{v}_{{0}}\frac{1-{{\text{e}}^{-\gamma t}}}{\gamma}\right]=0 \end{array}\end{eqnarray} \tag{ 25 }$$

and the conditional second moment of the displacement is

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {} & \,\mathbb{E}\left(|\boldsymbol{x}(t)-\boldsymbol{x}_{{0}}{{|}^{2}}\,|\,\boldsymbol{x}_{{0}},\boldsymbol{v}_{{0}}\right)={\int}^{}|\boldsymbol{x}-\boldsymbol{x}_{{0}}{{|}^{2}}p\,\left(\boldsymbol{x},t\,|\,\boldsymbol{x}_{{0}},\boldsymbol{v}_{{0}}\right)\,\text{d}\boldsymbol{x} \\ {} & = & \,\frac{|\boldsymbol{v}_{{0}}{{|}^{2}}}{{{\gamma}^{2}}}{{\left(1-{{\text{e}}^{-\gamma t}}\right)}^{2}}+\frac{3kT}{m{{\gamma}^{2}}}\left(2\gamma t-3+4{{\text{e}}^{-\gamma t}}-{{\text{e}}^{-2\gamma t}}\right), \end{array}\end{eqnarray} \tag{ 26 }$$

which is independent of $\boldsymbol{x}_{{0}}$. Using the Maxwell distribution of velocities (14), we find that the unconditional second moment is

$$\begin{eqnarray}\begin{array}{*{35}{l}} \mathbb{E}|\boldsymbol{x}(t)-\boldsymbol{x}_{{0}}{{|}^{2}}\,\,& = \,{{\mathbb{E}}_{\boldsymbol{x}_{{0}}}}{{\mathbb{E}}_{\boldsymbol{v}_{{0}}}}\left(|\boldsymbol{x}(t)-\boldsymbol{x}_{{0}}{{|}^{2}}\,|\,\boldsymbol{x}_{{0}},\boldsymbol{v}_{{0}}\right) \\ & = \frac{3kT}{m{{\gamma}^{2}}}{{\left(1-{{\text{e}}^{-\gamma t}}\right)}^{2}}\,\,+\,\frac{3kT}{m{{\gamma}^{2}}}\left(2\gamma t-3+4{{\text{e}}^{-\gamma t}}-{{\text{e}}^{-2\gamma t}}\right) \\ & = \frac{6kT}{m{{\gamma}^{2}}}\left(\gamma t-1+{{\text{e}}^{-\gamma t}}\right). \end{array}\end{eqnarray} \tag{ 27 }$$

The long time asymptotics of $\mathbb{E}|\boldsymbol{x}(t)-\boldsymbol{x}_{{0}}{{|}^{2}}$ is found from (27) to be

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathbb{E}|\boldsymbol{x}(t)-\boldsymbol{x}_{{0}}{{|}^{2}}\sim \frac{6kT}{m\gamma}t=\frac{kT}{ma\eta}t~\text{for}~t\gamma \gg 1; \end{array}\end{eqnarray} \tag{ 28 }$$

that is, the displacement variance of each component is asymptotically $3kT/ma\eta$. It was this fact that was verified experimentally by Perrin [21]. The one-dimensional diffusion coefficient, as defined in (4), is therefore given by $D=kT/6ma\eta$.

Equation (27) implies that the short time asymptotics of $\mathbb{E}|\boldsymbol{x}(t)-\boldsymbol{x}_{{0}}{{|}^{2}}$ is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathbb{E}|\boldsymbol{x}(t)-\boldsymbol{x}_{{0}}{{|}^{2}}\sim \frac{3kT}{m}{{t}^{2}}=\langle |\boldsymbol{v}_{{0}}{{|}^{2}}\rangle {{t}^{2}}. \end{array}\end{eqnarray} \tag{ 29 }$$

This result was first obtained by Smoluchowski.

To reconcile the Einstein and the Langevin approaches, we have to show that for two disjoint time intervals, $\left({{t}_{1}},{{t}_{2}}\right)$ and $\left({{t}_{3}},{{t}_{4}}\right)$, in the limit $\gamma \to \infty$, the increments ${{ \Delta }_{1}}\boldsymbol{x}=\boldsymbol{x}\left({{t}_{2}}\right)-\boldsymbol{x}\left({{t}_{1}}\right)$ and ${{ \Delta }_{3}}\boldsymbol{x}=\boldsymbol{x}\left({{t}_{4}}\right)-\boldsymbol{x}\left({{t}_{3}}\right)$ are independent zero mean Gaussian variables with variances proportional to the time increments. Equation (23) implies that in the limit $\gamma \to \infty$ the increments ${{ \Delta }_{1}}\boldsymbol{x}$ and ${{ \Delta }_{3}}\boldsymbol{x}$ are zero mean Gaussian variables and (28) shows that the variance of an increment is proportional to the time increment.

To show that the increments are independent, we use (20) in (23) to obtain

$$\begin{eqnarray}\begin{array}{*{35}{l}} \underset{\gamma \to \infty}{{\lim}}\,{{\gamma}^{2}}\langle {{ \Delta }_{1}}\boldsymbol{x}\centerdot {{ \Delta }_{3}}\boldsymbol{x}\rangle & = \,\underset{\gamma \to \infty}{{\lim}}\,{\int}_{{{t}_{1}}}^{{{t}_{2}}}{\int}_{{{t}_{3}}}^{{{t}_{4}}}\langle \gamma \boldsymbol{v}\left({{s}_{1}}\right)\centerdot \gamma \boldsymbol{v}\left({{s}_{2}}\right)\rangle \,\text{d}{{s}_{1}}\,\text{d}{{s}_{2}} \\ & = \,\frac{2q}{{{m}^{2}}}{\int}_{{{t}_{1}}}^{{{t}_{2}}}{\int}_{{{t}_{3}}}^{{{t}_{4}}}\delta \left({{s}_{2}}-{{s}_{1}}\right)\,\text{d}{{s}_{1}}\,\text{d}{{s}_{2}}=0. \end{array}\end{eqnarray} \tag{ 30 }$$

As is well-known [30], uncorrelated Gaussian variables are independent. This reconciles the Einstein and Langevin theories of Brownian motion in liquid.

Introducing the dimensionless variables $s=\gamma t$ and $\boldsymbol{\xi }(s)=\sqrt{m/6kT}\gamma \boldsymbol{x}(t)$, we find from (27) that

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\gamma \to \infty}{{\lim}}\,\mathbb{E}|\boldsymbol{\xi }(s)-\boldsymbol{\xi }(0){{|}^{2}}\,\,=s-1+{{\text{e}}^{-s}}\sim s\,\text{for}\,s\gg 1 \end{array}\end{eqnarray} \tag{ 31 }$$

and from (30) that

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{\gamma \to \infty}{{\lim}}\,\mathbb{E}\left[{{ \Delta }_{1}}\boldsymbol{\xi }\centerdot {{ \Delta }_{3}}\boldsymbol{\xi }\right]=0. \end{array}\end{eqnarray} \tag{ 32 }$$

Equations (31) and (32) explain (in the context of Langevin's description) Einstein's quoted assumption that ' ... the movement of one and the same particle after different intervals of time (are) mutually independent processes, so long as we think of these intervals of time as being chosen not too small'.

Starting with the Smoluchowski equation in dimension $d\geqslant 1$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \overset{\centerdot}{{\boldsymbol{x}}}\,(t)=\boldsymbol{a}\left(\boldsymbol{x}(t),t\right)+\boldsymbol{b}\left(\boldsymbol{x}(t),t\right)\,\overset{\centerdot}{{\boldsymbol{w}}}\,(t), \end{array}\end{eqnarray} \tag{ 38 }$$

where $\boldsymbol{x}$ and $\boldsymbol{a}\left(\boldsymbol{x},t\right)$ are d-dimensional vectors, $\boldsymbol{w}(t)$ is an m-dimensional Brownian motion, and $\boldsymbol{b}\left(\boldsymbol{x},t\right)$ is a $d\times m$ matrix, approximate trajectories of (38) can be constructed on any time interval s  <  t  <  T by the Euler scheme

$$\begin{eqnarray}\begin{array}{*{35}{l}} \boldsymbol{x}\left(t+ \Delta t\right)= & \,\boldsymbol{x}(t)+\boldsymbol{a}\left(\boldsymbol{x}(t),t\right) \Delta t+\boldsymbol{b}\left(\boldsymbol{x}(t),t\right)\, \Delta \boldsymbol{w}\left(t, \Delta t\right) \end{array}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \boldsymbol{x}_{{N}}(s)= & \boldsymbol{x}_{{0}}. & \end{array}\end{eqnarray}$$

Boundary behavior can be imposed on the trajectories of (39) in a given domain D. For example, all trajectories of (39) can be instantaneously terminated when they exit D. In this case the boundary $\partial D$ is called an absorbing boundary. The trajectories can be instantaneously reflected at $\partial D$ back into D according to a given reflection law; in this case $\partial D$ is called a reflecting boundary. They can also be either instantaneously terminated with a given probability or instantaneously reflected; then $\partial D$ is called a partially reflecting boundary.

It can be shown that the trajectories of (39) with a given boundary behavior converge to a limit as $\Delta t\to 0$. The limit is defined as the solution of the Smoluchowski equation (38) with the given boundary behavior. (See [18] and [55] for details and for other types of boundary behavior.)

We assume that in the one-dimensional case $b(x,t)>\delta >0$ for some constant δ. We assume for now that a(x,t) and b(x,t) are deterministic functions. To construct the pdf of a trajectory of (39), we note that the pdf of x_N(t) can be expressed explicitly for t on the lattice, because (39), written as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta w(t)=\frac{{{x}_{N}}\left(t+ \Delta t\right)-{{x}_{N}}(t)-a\left({{x}_{N}}(t),t\right)}{b\left({{x}_{N}}(t),t\right)}, \end{array}\end{eqnarray} \tag{ 40 }$$

means that for all t on the lattice the expressions on the right-hand side of (40) are independent, identically distributed Gaussian variables. It follows that the pdf of the entire Euler trajectory is the product

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \,p\left({{x}_{1}},{{t}_{1}};{{x}_{2}},{{t}_{2}};\ldots ;{{x}_{n}},{{t}_{n}}\right) \\ {} & =\,\underset{k=1}{\overset{n}{\prod}}\,{{\left[2\pi {{b}^{2}}\left({{x}_{k-1}},{{t}_{k-1}}\right) \Delta t\right]}^{-1/2}}\exp \left\{-\frac{{{\left[{{x}_{k}}-{{x}_{k-1}}-a\left({{x}_{k-1}},{{t}_{k-1}}\right) \Delta t\right]}^{2}}}{2{{b}^{2}}\left({{x}_{k-1}},{{t}_{k-1}}\right) \Delta t}\right\}. \end{array}\end{eqnarray} \tag{ 41 }$$

Setting x_n  =  x and integrating over $\mathbb{R}$ with respect to all intermediate points ${{x}_{1}},{{x}_{2}},\,\ldots,{{x}_{n-1}}$, we find from (41) that the transition pdf of the trajectory satisfies on the lattice the recurrence relation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{N}}\left(x,t+ \Delta t\,|\,{{x}_{0}}\right)={\int}_{\mathbb{R}}\frac{{{p}_{N}}\left(\,y,t\,|\,{{x}_{0}}\right)\,\text{d}y}{\sqrt{2\pi \Delta t}\,b\left(\,y,t\right)}\exp \left\{-\frac{{{\left[x-y-a\left(\,y,t\right) \Delta t\right]}^{2}}}{2{{b}^{2}}\left(\,y,t\right) \Delta t}\right\}. \end{array}\end{eqnarray} \tag{ 42 }$$

The solution of the integral equation (42) is called Wiener's discrete path integral. Its limit as $N\to \infty$ is called Wiener's path integral.

For $\Delta t=(t-s)/N$, in the limit $N\to \infty$, the pdf ${{p}_{N}}\left(x,t\,|\,{{x}_{0}}\right)$ of the solution ${{x}_{N}}\left(t,\mathfrak{w}\right)$ of (39) converges to the solution $p\,\left(x,t\,|\,{{x}_{0}}\right)$ of (42), where $\mathfrak{w}$ is the entire discrete path of the Brownian motion w(t) on the lattice. Expansion of ${{p}_{N}}\left(x,t\,|\,{{x}_{0}}\right)$ in (42) shows that the pdf $p\left(x,t\,|\,{{x}_{0}}\right)$ is also the solution of the initial value problem

$$\begin{eqnarray}&&\frac{\partial p\,\left(\,y,t\,|\,x,s\right)}{\partial t}=\frac{1}{2}\frac{{{\partial}^{2}}\left[{{b}^{2}}\left(\,y,t\right)p\,\left(\,y,t\,|\,x,s\right)\right]}{\partial {{y}^{2}}}-\frac{\partial \left[a\left(\,y,t\right)p\,\left(\,y,t\,|\,x,s\right)\right]}{\partial y}\end{eqnarray} \tag{ 43 }$$

with the initial condition

$$\begin{eqnarray}&&\underset{t\downarrow s}{{\lim}}\,p\,\left(\,y,t\,|\,x,s\right)=\delta \left(\,y-x\right).\end{eqnarray} \tag{ 44 }$$

Equation (43) is called the Fokker–Planck–(Smoluchowski) equation (FPE) (see [55]).

In the Smoluchowski equation (38) for $d=n\geqslant 1$, the $n\times m$ matrix

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{B}\left(\boldsymbol{x},t\right)={{\left\{{{b}^{ij}}\left(\boldsymbol{x},t\right)\right\}}_{n\times m}} \end{array}\end{eqnarray} \tag{ 45 }$$

is called the noise matrix, and $\boldsymbol{\sigma }\left(\boldsymbol{x},t\right)=\frac{1}{2}\boldsymbol{B}\left(\boldsymbol{x},t\right)\boldsymbol{B}^{{T}}\left(\boldsymbol{x},t\right)$ is called the diffusion matrix. The operator

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{L}}_{\boldsymbol{y}}}p=\underset{i=1}{\overset{d}{\sum}}\,\frac{\partial}{\partial {{y}^{i}}}\left\{\underset{j=1}{\overset{d}{\sum}}\,\frac{\partial}{\partial {{y}^{j}}}{{\sigma}^{ij}}\left(\boldsymbol{y},t\right)p-{{a}^{i}}\left(\boldsymbol{y},t\right)p\right\}, \end{array}\end{eqnarray} \tag{ 46 }$$

is called the Fokker–Planck operator, or the forward Kolmogorov operator.

As in the one-dimensional case, the pdf $p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)$ that $\boldsymbol{x}(t)=\boldsymbol{y}$, given that $\boldsymbol{x}(s)=\boldsymbol{x}$, is the solution of the Fokker–Planck–Smoluchowski initial value problem

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{\partial p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)}{\partial t}= & \,{{\mathcal{L}}_{\boldsymbol{y}}}p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)~\text{for}~\boldsymbol{x},\boldsymbol{y}\in {{\mathbb{R}}^{n}},~t>s, \end{array}\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \underset{t\downarrow s}{{\lim}}\,p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)= & \,\delta \left(\boldsymbol{x}-\boldsymbol{y}\right). \end{array}\end{eqnarray} \tag{ 48 }$$

It can be written as the conservation law by defining the probability flux density vector

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{J}^{i}}\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)=-\underset{j=1}{\overset{d}{\sum}}\,\frac{\partial}{\partial {{y}^{j}}}{{\sigma}^{ij}}\left(\boldsymbol{y},t\right)p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)-{{a}^{i}}\left(\boldsymbol{y},t\right)p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right) \end{array}\end{eqnarray} \tag{ 49 }$$

and writing (47) in the divergence form

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\partial p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)}{\partial t}=-\text{di}{{\text{v}}_{\boldsymbol{y}}}\,\boldsymbol{J}\left(\boldsymbol{y},t\,|\boldsymbol{x},s\right). \end{array}\end{eqnarray} \tag{ 50 }$$

The probability density at time t and at point $\boldsymbol{x}$ can be represented for any domain $\Omega$ by the limit as $N\to \infty$ of

$$\begin{eqnarray}&&\begin{array}{l} {\rm Pr}\left\{{\boldsymbol{x}}_N(t_{1,N})\in\Omega,{\boldsymbol{x}}_N(t_{2,N})\in\Omega,\dots, {\boldsymbol{x}}_N(t)=\boldsymbol{x}, t\leq T\leq t+\Delta t\,\vert \,\boldsymbol{x}(0)=\boldsymbol{y}\right\}\\ ={\left[\right.}\int_{\Omega} \int_{\Omega}\cdots \int_{\Omega}\,\prod_{j=1}^{N} \frac{{\rm d}{\boldsymbol{y}}_j}{\sqrt{(2\pi \Delta t)^n\det\boldsymbol{\sigma}(\boldsymbol{x})(t_{j-1,N}))}} \\ \qquad \times \exp \left\{ -\frac{1}{2\Delta t} \left[\boldsymbol{\boldsymbol{y}}_j-\boldsymbol{x}(t_{j-1,N})- \boldsymbol{a}({\boldsymbol{x}}(t_{j-1,N}))\Delta t \right]^T\boldsymbol{\sigma}^{-1}(\boldsymbol{x}(t_{j-1,N}))\right. \\ \left. \qquad \times\left[{\boldsymbol{y}}_j-\boldsymbol{x}(t_{j-1,N})-\boldsymbol{a}(\boldsymbol{x}(t_{j-1,N}))\Delta t \right]\right\}, \end{array}\end{eqnarray} \tag{ 51 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Delta t=\frac{t}{N},\quad {{t}_{j,N}}=j \Delta t, \end{array}\end{eqnarray}$$

and

$$\begin{eqnarray}&&\boldsymbol{x}\left({{t}_{0,N}}\right)=\boldsymbol{y}\end{eqnarray}$$

in the product. The limit is the Wiener integral defined by the stochastic differential equation (38) with appropriate boundary condition. In the limit $N\to \infty$ the integral (51) converges to the solution of the Fokker–Planck equation (50) in $\Omega$. Similarly, expanding the path integral with respect to $\boldsymbol{x}$ and s, we find that $p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)$ satisfies the backward Kolmogorov equation [18]

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{L}_{\boldsymbol{x}}^{\ast}p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)=\underset{i=1}{\overset{n}{\sum}}\,\underset{j=1}{\overset{n}{\sum}}\,{{\sigma}^{ij}}\left(\boldsymbol{x},s\right)\frac{{{\partial}^{2}}p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)}{\partial {{x}^{i}}\partial {{x}^{j}}}+\underset{i=1}{\overset{n}{\sum}}\,{{a}^{i}}\left(\boldsymbol{x},s\right)\frac{\partial p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)}{\partial {{x}^{i}}} \end{array}\end{eqnarray} \tag{ 52 }$$

with the terminal condition

$$\begin{eqnarray}&&\underset{t\uparrow s}{{\lim}}\,p\,\left(\,\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)=\delta \left(\boldsymbol{y}-\boldsymbol{x}\right).\end{eqnarray} \tag{ 53 }$$

The Fokker–Planck equation can be generalized when a killing measure is added to the dynamics. A killing measure represents the probability per unit time and unit length to terminate a trajectory at a given point at a given time. Mixed boundary conditions and Fokker–Planck equation for the survival probability with killing is discussed in [14].

A d-dimensional Markov process $\boldsymbol{x}(t)$ is called a diffusion process with (deterministic) drift vector field $\boldsymbol{a}\left(\boldsymbol{x},t\right)$ and (deterministic) diffusion matrix $\boldsymbol{\sigma }\left(\boldsymbol{x},t\right)$, if it has continuous trajectories,

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \underset{ \Delta t\to 0}{{\lim}}\,\frac{1}{ \Delta t}\mathbb{E}\left\{\boldsymbol{x}\left(t+ \Delta t\right)-\boldsymbol{x}(t)\,|\,\boldsymbol{x}(t)=\boldsymbol{x}\right\}=\boldsymbol{a}\left(\boldsymbol{x},t\right) \end{array}\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&\underset{ \Delta t\to 0}{{\lim}}\,\frac{1}{ \Delta t}\mathbb{E}\left\{\left[{{x}^{i}}\left(t+ \Delta t\right)-{{x}^{i}}(t)\right]\left[{{x}^{j}}\left(t+ \Delta t\right)-{{x}^{j}}(t)\right]\,|\,\boldsymbol{x}(t)=\boldsymbol{x}\right\}\,={{\sigma}^{ij}}\left(\boldsymbol{x},t\right)\end{eqnarray} \tag{ 55 }$$

for $i,j=1,2,\ldots,d$, and for some $\delta >0$

$$\begin{eqnarray}&&\underset{ \Delta t\to 0}{{\lim}}\,\frac{1}{ \Delta t}\mathbb{E}\left\{|\boldsymbol{x}\left(t+ \Delta t\right)-\boldsymbol{x}(t)|{{\,}^{2+\delta}}\,|\,\,\boldsymbol{x}(t)=\boldsymbol{x}\right\}=0.\end{eqnarray}$$

It can be shown [18] that under mild regularity assumptions on the coefficient the solution of the stochastic differential equation (SDE)

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \overset{\centerdot}{{\boldsymbol{x}}}\,(t)=\boldsymbol{a}\left(\boldsymbol{x}(t),t\right)+\boldsymbol{B}\left(\boldsymbol{x}(t),t\right)\,\overset{\centerdot}{{\boldsymbol{w}}}\,(t) \end{array}\end{eqnarray} \tag{ 56 }$$

is a diffusion process with drift field $\boldsymbol{a}\left(\boldsymbol{x},t\right)$ and diffusion matrix $\boldsymbol{\sigma }\left(\boldsymbol{x},t\right)=\frac{1}{2}\boldsymbol{B}\left(\boldsymbol{x},t\right)\boldsymbol{B}^{{T}}\left(\boldsymbol{x},t\right)$. Also a partial converse is true: given sufficiently nice drift field $\boldsymbol{a}\left(\boldsymbol{x},t\right)$ and a strictly positive definite diffusion matrix $\boldsymbol{\sigma }\left(\boldsymbol{x},t\right)$ of a diffusion process, there is a noise matrix $\boldsymbol{B}\left(\boldsymbol{x},t\right)$ and a Brownian motion $\boldsymbol{w}(t)$ such the process is a solution of the SDE (56) [18, 19].

Equations (54) and (55) reconstruct the SDE from its trajectories. The discrete version of these equations can be used as approximations to the coefficients, when a sufficiently large set of trajectory fragments is given (see [56, 58, 59]). Indeed, a large number (tens of thousands) of short single particle trajectories (SPTs), collected by super-resolution methods, at tens of nanometer precision for motion occurring in cell are ideal data to recover drift and diffusion tensors [67]. But the sampled molecular process cannot be directly modeled by the overdamped Langevin equation, which describes diffusion on the microscopic level, because an additional noise of localization is added in tracking the motion (algorithm and point spread function error). Using a stochastic model of the acquired data to calibrate the model, it is possible to distinguish the perturbation added to the physical motion (stochastic equation) [58]. The denoising procedure of the SPTs reveals that effective diffusion coefficient contains the divergence of deterministic drift component and also provides a criteria to differentiate trapped stochatic particle from immobile ones [58].

If the Smoluchowski trajectories are terminated at the boundary $\partial D$ of a given domain, the pdf vanishes on $\partial D$. Thus equations (47) and (48) have to be supplemented by the boundary condition

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)=0~\text{for}~\boldsymbol{y}\in \partial D,~\boldsymbol{x}\in D,~t>s. \end{array}\end{eqnarray} \tag{ 57 }$$

In this case, the first passage time to the boundary $\partial D$ is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\tau}_{D}}=\inf \left\{t:\,\boldsymbol{x}(t)\notin D\right\}. \end{array}\end{eqnarray} \tag{ 58 }$$

Thus, the probability that the Smoluchowski trajectory is still in D at time t  >  s, given that at time s  <  t it started at $\boldsymbol{x}\in D$, is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{Pr}\left\{\boldsymbol{x}(t)\in D\,|\,\boldsymbol{x}(s)=\boldsymbol{x}\right\}={{{\int}^{}}_{D}}p\left(\,\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)\,\text{d}\boldsymbol{y} \end{array}.\end{eqnarray} \tag{ 59 }$$

Actually, $\text{Pr}\left\{\boldsymbol{x}(t)\in D\,|\,\boldsymbol{x}(s)=\boldsymbol{x}\right\}=\text{Pr}\left\{{{\tau}_{D}}>t\,|\,\boldsymbol{x}(s)=\boldsymbol{x}\right\},$ which is the survival probability at time t  >  s of the Smoluchowski trajectory that started at $\boldsymbol{x}$ at time s  <  t.

If the Smoluchowski trajectories are reflected at the boundary $\partial D$ of a given domain such that the normal flux at the boundary vanishes [16], then the boundary condition becomes

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{J}\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)\centerdot \boldsymbol{\nu }\left(\boldsymbol{y}\right)=0\,\text{for}\,\,\boldsymbol{y}\in \partial D,\boldsymbol{x}\in D,\,t>s, \end{array}\end{eqnarray} \tag{ 60 }$$

where $\boldsymbol{\nu }\left(\boldsymbol{y}\right)$ is the unit outer normal vector at $\boldsymbol{y}\in \partial D$.

When the boundary $\partial D$ absorbs the Smoluchowski trajectories, the expected MFPT to the boundary of Smoluchowski trajectories in D is found by integrating the survival probability to obtain the MFPT as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathbb{E}\left[{{\tau}_{D}}\,|\,\boldsymbol{x}(s)=\boldsymbol{x}\right]={\int}_{s}^{\infty}{\int}_{D}p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)\,\text{d}\boldsymbol{y}\,\text{d}t. \end{array}\end{eqnarray} \tag{ 61 }$$

Setting $u\left(\boldsymbol{x},s\right)=\mathbb{E}\left[{{\tau}_{D}}\,|\,\boldsymbol{x}(s)=\boldsymbol{x}\right]$, using the backward Kolmogorov equation and Green's identity, we obtain the backward boundary value problem

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{L}}^{\ast}}_{\boldsymbol{x}}u\left(\boldsymbol{x},s\right)+\frac{\partial u\left(\boldsymbol{x},s\right)}{\partial s}= & -1~\text{for}~\boldsymbol{x}\in D, \end{array}\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} u\left(\boldsymbol{x},s\right)= & \,0~\text{for}~\boldsymbol{x}\in \partial D. \end{array}\end{eqnarray} \tag{ 63 }$$

Note that if the coefficients $\boldsymbol{\sigma }$ and $\boldsymbol{a}$ are time-independent, then (52), (63) reduce to the time-homogeneous elliptic Pontryagin–Andronov–Vitt (PAV) [17] boundary value problem in D

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{L}}^{\ast}}_{\boldsymbol{x}}u\left(\boldsymbol{x}\right)= & -1~\text{for}~\boldsymbol{x}\in D, \end{array}\end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} u\left(\boldsymbol{x}\right)= & \,0~\text{for}~\boldsymbol{x}\in \partial D. \end{array}\end{eqnarray} \tag{ 65 }$$

If the boundary is reflecting, then the absorbing boundary condition (65) for the PAV equation (64) is changed to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{J}\left(\boldsymbol{x}\right)\centerdot \boldsymbol{\nu }\left(\boldsymbol{x}\right)=\,0~\text{for}~\boldsymbol{x}\in \partial D. \end{array}\end{eqnarray} \tag{ 66 }$$

If $\partial D$ is absorbing on a part $\partial {{D}_{a}}$ and reflecting on the remaining part $\partial {{D}_{r}}-\partial D-\partial {{D}_{a}}$, then the boundary conditions for the PAV equation is absorbing on $\partial {{D}_{a}}$ and reflecting on $\partial {{D}_{b}}$.

Note further that the MFPT to the boundary is given by (61), where $p\left(\boldsymbol{y},t\,|\,\boldsymbol{x},s\right)$ is the solution of the FPE with absorbing boundary conditions also in the case of reflecting boundaries, because prior to reaching the boundary the Smoluchowski trajectories are independent of boundary behavior.

We summarize here the asymptotic formulas of the MFPT (67) when the domain $\Omega$ is in the plane and absorbing boundary is a small sub arc $\partial {{ \Omega }_{a}}$ (of length a) of the boundary $\partial \Omega$. We have reviewed the mathematical method and the analysis in [60, 62]. There is little intuition behind these formulas and it is not fruitful to guess what there are. It is indeed hard to tell in advance how the geometry enters into the formulas. The recipe we adopted is to follow the analytical derivations that reveal how local and global structures, smoothness or not, local curvature controls the narrow escape time.

• 1.  
  When $\partial {{ \Omega }_{a}}$ is a sub-arc of a smooth boundary, the MFPT from any point $\boldsymbol{x}$ in $\Omega$ to $\partial {{ \Omega }_{a}}$ is denoted ${{\bar{\tau}}_{\boldsymbol{x}\to \partial {{ \Omega }_{a}}}}$. For

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \varepsilon =\frac{\pi |\partial {{ \Omega }_{a}}|}{|\partial \Omega |}=\frac{\pi a}{|\partial \Omega |}\ll 1 \end{array}\end{eqnarray} \tag{ 77 }$$

  the MFPT is independent of $\boldsymbol{x}$ outside a small vicinity of $\partial {{ \Omega }_{a}}$ (called a boundary layer). Thus for $\boldsymbol{x}\in \Omega$, outside a boundary layer near $\partial {{ \Omega }_{a}}$,

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\bar{\tau}}}_{\boldsymbol{x}\to \partial {{ \Omega }_{a}}}}=\frac{| \Omega |}{\pi D}\ln \frac{1}{\varepsilon}+O(1), \end{array}\end{eqnarray} \tag{ 78 }$$

  where O(1) depends on the initial distribution of $\boldsymbol{x}$ [35–57]. This result was derived independently using matched asymptotic technics and Green's function method.If $\Omega$ is a disc of radius R, then for $\boldsymbol{x}$ at the center of the disk (figure 2(A)),

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\bar{\tau}}}_{\boldsymbol{x}\to \partial {{ \Omega }_{a}}}}=\frac{{{R}^{2}}}{D}\left[\log \frac{R}{a}+2\log 2+\frac{1}{4}+O\left(\varepsilon \right)\right], \end{array}\end{eqnarray}$$

  and averaging with respect to a uniform distribution of $\boldsymbol{x}$ in the disk [14]

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}=\frac{{{R}^{2}}}{D}\left[\log \frac{R}{a}+2\log 2+\frac{1}{8}+O\left(\varepsilon \right)\right]. \end{array}\end{eqnarray}$$

  This result was obtained from generalizing Sneddon's method for mixed boundary value problem. The method is based on Abel's transformation [40–42]. The flux through a hole in a smooth wall on a flat membrane surface is regulated by the area $| \Omega |$ inside the wall, the diffusion coefficient D, and the aspect ratio ε (77). In the case of Brownian motion on a sphere of radius R the MFPT to an absorbing circle centered on the north-south axis near the south pole with small radius $a=R\sin \delta /2$ is given by

  $$\begin{eqnarray}&&\bar{\tau}=\frac{2{{R}^{2}}}{D}\log \frac{\sin \frac{\theta}{2}}{\sin \frac{\delta}{2}},\end{eqnarray} \tag{ 79 }$$

  where θ is the angle between $\boldsymbol{x}$ and the south-north axis of the sphere (figure 2(B)).
• 2.  
  If the absorbing window is located at a corner of angle α, then

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}=\frac{| \Omega {{|}_{g}}}{D\alpha}\left[\log \frac{1}{\varepsilon}+O(1)\right], \end{array}\end{eqnarray} \tag{ 80 }$$

  where $| \Omega {{|}_{g}}$ is the surface area of the domain on the curved surface, calculated according to the Riemannian metric on the surface [40]. Formula (80) indicates that control of flux is regulated also by the access to the absorbing window afforded by the angle of the corner leading to the window (figure 2(C)). This formula was obtained using a conformal map sending a corner to a flat line.
• 3.  
  If the absorbing window is located at a cusp, then $\bar{\tau}$ grows algebraically, rather than logarithmically. Thus, in the domain bounded between two tangent circles, the expected lifetime is

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}=\frac{| \Omega |}{\left({{d}^{-1}}-1\right)D}\left(\frac{1}{\varepsilon}+O(1)\right), \end{array}\end{eqnarray} \tag{ 81 }$$

  where d  <  1 is the ratio of the radii [42] (figure 2(F)). Formula (81) indicates that a drastic reduction of flux can be achieved by putting an obstacle that limits the access to the absorbing window by forming a cusp-like passage. This formula was derived using the exponential conformal map.
• 4.  
  When $\partial {{ \Omega }_{a}}$ (of length a) is located at the end of a narrow neck with radius of curvature R_c, the MFPT is given in [14, 63] as (figures 2(G) and (I))

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}=\frac{| \Omega |}{4D\sqrt{2a/{{R}_{c}}}}\left(1+O(1)\right)~\text{for}~a\ll |\partial \Omega |. \end{array}\end{eqnarray} \tag{ 82 }$$

  This formula is derived by a new method that uses a Mobius transformation to resolve the cusp singularity [14, 63]. The boundary layer at the cusp is sent to a banana shaped domain. Asymptotic formula for a general cusp with an arbitrary power law are not known.For a surface of revolution generated by rotating the curve about its axis of symmetry [63], we use the representation of the generating curve

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} y=r(x),~ \Lambda <x<0~ \end{array}\end{eqnarray}$$

  where the x-axis is horizontal with $x= \Lambda$ at the absorbing end $\boldsymbol{A}\boldsymbol{B}$. We assume that the parts of the curve that generate the funnel have the form

  $$\begin{eqnarray}\begin{array}{*{35}{l}} r(x) & =O\left(\sqrt{|x|}\right)\,\text{near}\,x=0 \\ r(x) & =a+\frac{{{\left(x- \Lambda \right)}^{1+\nu}}}{\nu \left(1+\nu \right){{\ell}^{\nu}}}\left(1+o(1)\right)\,\text{for}\,\nu >0\,\text{near}\,x= \Lambda, \end{array}\end{eqnarray} \tag{ 83 }$$

  where $a=\frac{1}{2}\overline{\boldsymbol{A}\boldsymbol{B}}=\varepsilon /2$ is the radius of the gap, and the constant $\ell$ has dimension of length. For $\nu =1$ the parameter $\ell$ is the radius of curvature R_c at $x= \Lambda$. The MFPT from the head to the absorbing end $\boldsymbol{A}\boldsymbol{B}$ is given by

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}\sim \frac{\mathcal{S}\left(\Lambda \right)}{2D}\frac{{{\left(\frac{\ell}{\left(1+\nu \right)a}\right)}^{\nu /1+\nu}}{{\nu}^{1/1+\nu}}}{\sin \frac{\nu \pi}{1+\nu}}, \end{array}\end{eqnarray} \tag{ 84 }$$

  where $\mathcal{S}$ is the entire unscaled area of the surface. In particular, for $\nu =1$ the MFPT (84) reduces to

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}\sim \frac{\mathcal{S}}{4D\sqrt{a/2\ell}}. \end{array}\end{eqnarray} \tag{ 85 }$$
• 5.  
  When a bulky head is connected to an essentially one-dimensional strip (or cylinder) of small radius a and length L, as is the case of a neuronal spine membrane (figure 2(D)). The connection of the head to the neck can be at an angle or by a smooth funnel. The boundary of the domain reflects Brownian trajectories and only the end of the cylinder $\partial {{ \Omega }_{a}}$ absorbs them. The domain ${{ \Omega }_{1}}$ is connected to the cylinder at an interface $\partial {{ \Omega }_{i}}$, which in this case is an interval AB. The MFPT from $\boldsymbol{x}\in {{ \Omega }_{1}}$ to $\partial {{ \Omega }_{a}}$ is given by

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\bar{\tau}}}_{\boldsymbol{x}\to \partial {{ \Omega }_{a}}}}={{{\bar{\tau}}}_{\boldsymbol{x}\to \partial {{ \Omega }_{i}}}}+\frac{{{L}^{2}}}{2D}+\frac{|{{ \Omega }_{1}}|L}{|\partial {{ \Omega }_{a}}|D}. \end{array}\end{eqnarray} \tag{ 86 }$$

  The flux dependence on the neck length is quite strong. This formula is derived using the additive property of the MFPT [66].
• 6.  
  A dumbbell-shaped domain (of type (VI)) consists of two compartments ${{ \Omega }_{1}}$ and ${{ \Omega }_{3}}$ and a connecting neck ${{ \Omega }_{2}}$ that is effectively one-dimensional (figure 2(J)), or in a similar domain with a long neck. A Brownian trajectory that hits the segment $\boldsymbol{A}\boldsymbol{B}$ in the center of the neck ${{ \Omega }_{2}}$ is equally likely to reach either compartment before the other; thus $\boldsymbol{A}\boldsymbol{B}$ is the stochastic separatrix (SS). Therefore the mean time to traverse the neck from compartment ${{ \Omega }_{1}}$ to compartment ${{ \Omega }_{3}}$ is asymptotically twice the MFPT ${{\bar{\tau}}_{{{ \Omega }_{1}}\to SS}}$. Neglecting, as we may, the mean residence time of a Brownian trajectory in ${{ \Omega }_{2}}$ relative to that in ${{ \Omega }_{1}}$ or in ${{ \Omega }_{3}}$ we can write the transition rates from ${{ \Omega }_{1}}$ to the ${{ \Omega }_{3}}$ and vv as

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\lambda}_{{{ \Omega }_{1}}\to {{ \Omega }_{3}}}}=\frac{1}{2{{{\bar{\tau}}}_{{{ \Omega }_{1}}\to SS}}},\quad {{\lambda}_{{{ \Omega }_{3}}\to {{ \Omega }_{1}}}}=\frac{1}{2{{{\bar{\tau}}}_{{{ \Omega }_{3}}\to SS}}}. \end{array}\end{eqnarray} \tag{ 87 }$$

  These rates can be found from explicit expressions for the flux into an absorbing window

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\lambda}_{1}}\sim \frac{1}{{\bar{\tau}}}, \end{array}\end{eqnarray} \tag{ 88 }$$

  where $\bar{\tau}$ is given in (86). Here ${{\bar{\tau}}_{\boldsymbol{x}\to \partial {{ \Omega }_{i}}}}$ is any one of the MFPTs given above, depending on the geometry of ${{ \Omega }_{1}}$ with L half the length of the neck and with $SS=\partial {{ \Omega }_{a}}$. The radii of curvature R_c,1 and R_c,3 at the two funnels may be different in ${{ \Omega }_{1}}$ and ${{ \Omega }_{3}}$. The smallest positive eigenvalue λ of the Neumann problem for the Laplace equation in the dumbbell is to leading order $\lambda =-\left({{\lambda}_{{{ \Omega }_{1}}\to {{ \Omega }_{3}}}}+{{\lambda}_{{{ \Omega }_{3}}\to {{ \Omega }_{1}}}}\right)$. For example, if the solid dumbbell consists of two general heads connected smoothly to the neck by funnels (see (93)), the two rates are given by

  $$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{1}{{{\lambda}_{{{ \Omega }_{1}}\to {{ \Omega }_{3}}}}}= & \sqrt{2}\left[{{\left(\frac{{{R}_{c,1}}}{a}\right)}^{3/2}}\frac{|{{ \Omega }_{1}}|}{{{R}_{c,1}}D}\right]\left(1+o(1)\right)+\frac{{{L}^{2}}}{4D}+\frac{|{{ \Omega }_{1}}|L}{\pi {{a}^{2}}D} \\ {} & {} \\ \frac{1}{{{\lambda}_{{{ \Omega }_{3}}\to {{ \Omega }_{1}}}}}= & \sqrt{2}\left[{{\left(\frac{{{R}_{c,3}}}{a}\right)}^{3/2}}\frac{|{{ \Omega }_{3}}|}{{{R}_{c,3}}D}\right]\left(1+o(1)\right)+\frac{{{L}^{2}}}{4D}+\frac{|{{ \Omega }_{3}}|L}{\pi {{a}^{2}}D} \end{array}\end{eqnarray} \tag{ 89 }$$

  (see [66]). Formulas (89) indicate that the unidirectional fluxes between the two compartments of a dumbbell-shaped domain can be controlled by the area (or surface area) of the two and by the type of obstacles to the access to the connecting neck. The equilibration rate in the dumbbell, λ, is thus controlled by the geometry.
• 7.  
  The mean time to escape through N well-separated absorbing windows of lengths a_j at the ends of funnels with radii of curvature ${{\ell}_{j}}$, respectively, in the boundary $\partial \Omega$ of a planar domain $\Omega$ is given by

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}=\frac{\pi | \Omega |}{2D\underset{j=1}{\overset{N}{\sum}}\,\sqrt{{{a}_{j}}/{{\ell}_{j}}}}\left(1+o(1)\right)~\text{for}~{{a}_{j}}/{{\ell}_{j}}\ll |\partial \Omega |. \end{array}\end{eqnarray} \tag{ 90 }$$

  The probability to escape through window i is given by

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{i}}=\frac{\sqrt{{{a}_{i}}/{{\ell}_{i}}}}{\underset{j=1}{\overset{N}{\sum}}\,\sqrt{{{a}_{j}}/{{\ell}_{j}}}}. \end{array}\end{eqnarray} \tag{ 91 }$$

  Formulas (90) and (91) are significant for diffusion in a network of compartments connected by narrow passages (e.g. on a membrane strewn with obstacles). The dependence of the MFPT $\bar{\tau}$ and of the transition probabilities p_i on the local geometrical properties of the compartments renders the effective diffusion tensor in the network position-dependent and can give rise to anisotropic diffusion.

Zoom In Zoom Out Reset image size

We now summarize the Narrow Escape Time formula in three-dimensions. The methods are the same as in two dimensions: Matched asymptotic or Greens function and conformal mapping to resolve cusp singularities. Indeed, the axial symmetry allows reducing the three- to two- dimensions and thus to use of conformal transformations [14].

• 1.  
  The MFPT to a circular absorbing window $\partial {{ \Omega }_{a}}$ of small radius a centered at $\mathbf{0}$ on the boundary $\partial \Omega$ is given by [43]

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\bar{\tau}}}_{\boldsymbol{x}\to \partial {{ \Omega }_{a}}}}=\frac{| \Omega |}{4aD\left[1+\frac{L\left(\mathbf{0}\right)+N\left(\mathbf{0}\right)}{2\pi}\,a\log a+o\left(a\log a\right)\right]}, \end{array}\end{eqnarray} \tag{ 92 }$$

  where $L\left(\mathbf{0}\right)$ and $N\left(\mathbf{0}\right)$ are the principal curvatures of the boundary at the center of $\partial {{ \Omega }_{a}}$. This formula was derived using the second order expansion of the Neuman-Green's function at the pole [43].
• 2.  
  The MFPT from the head of the solid of revolution, obtained by rotating the symmetric domain about its axis of symmetry, to a small absorbing window $\partial {{ \Omega }_{a}}$ at the end of a funnel (figure 2(H)) is given by

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}=\frac{1}{\sqrt{2}}{{\left(\frac{{{R}_{c}}}{a}\right)}^{3/2}}\frac{| \Omega |}{{{R}_{c}}D}\left(1+o(1)\right)~\text{for}~a\ll {{R}_{c}}, \end{array}\end{eqnarray} \tag{ 93 }$$

  where the R_c is the radius of curvature of the rotated curve at the end of the funnel [66].
• 3.  
  The MFPT from a point $\boldsymbol{x}$ in a bulky head $\Omega$ to an absorbing disk $\partial {{ \Omega }_{a}}$ of a small radius a at the end of a narrow neck of length L, connected to the head at an interface $\partial {{ \Omega }_{i}}$ is given by the connection formula (86). When the cylindrical neck is attached to the head at a right angle the interface $\partial {{ \Omega }_{i}}$ is a circular disk and ${{\bar{\tau}}_{\boldsymbol{x}\to \partial {{ \Omega }_{i}}}}$ is given by (92). When the neck is attached smoothly through a funnel, ${{\bar{\tau}}_{\boldsymbol{x}\to \partial {{ \Omega }_{i}}}}$ is given by (93).
• 4.  
  The mean time to escape through N well-separated absorbing circular windows or radii a_j at the ends of funnels with curvatures ${{\ell}_{j}}$, respectively, is given by

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{\tau}=\frac{1}{\sqrt{2}}\,\frac{| \Omega |}{D\underset{j=1}{\overset{N}{\sum}}\,{{\ell}_{j}}{{\left(\frac{{{a}_{j}}}{{{\ell}_{j}}}\right)}^{3/2}}}. \end{array}\end{eqnarray} \tag{ 94 }$$

  The exit probability through window i is given by

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{i}}=\frac{{{a}_{i}}^{3/2}\ell _{i}^{-1/2}}{\underset{j=1}{\overset{N}{\sum}}\,{{a}_{j}}^{3/2}\ell _{j}^{-1/2}}. \end{array}\end{eqnarray} \tag{ 95 }$$
• 5.  
  The principal eigenvalue of the Laplace equation in a dumbbell-shaped structure is given in item (vi), equations (87)–(89) above [66].
• 6.  
  The leakage flux through a circular hole of small radius a centered at $\mathbf{0}$ in the reflecting boundary is given by [43]

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{J}_{a}}=4aD{{u}_{0}}\left(\mathbf{0}\right)+O\left(\frac{{{a}^{2}}}{| \Omega {{|}^{2/3}}}\log \frac{a}{| \Omega {{|}^{1/3}}}\right), \end{array}\end{eqnarray} \tag{ 96 }$$

  where ${{u}_{0}}\left(\mathbf{0}\right)$ is the concentration of diffusers at the window in the same model without the absorbing window.

We refer the reader to the classical literature about the asymptotic formula for Narrow Escape time [61], the Dire Strait time (when escape occurs at a cusp boundary) [63] and the recent monograph [14] for applications in cellular biology.

The motion of monomer $\boldsymbol{R}_{{c}}$ of a Rouse polymer is related to the Fourier coefficients by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{R}_{{c}}=\underset{p=0}{\overset{N-1}{\sum}}\,\alpha _{p}^{c}\boldsymbol{u}_{{p}}, \end{array}\end{eqnarray} \tag{ 107 }$$

where $\alpha _{p}^{c}$ are described by relation (102) and $\boldsymbol{u}_{{p}}$ satisfy equations (105), which form an ensemble of Ornstein–Uhlenbeck processes, for which the variance is simply

$$\begin{eqnarray}\begin{array}{*{35}{l}} \sigma _{p}^{2}(t) & =\langle |\boldsymbol{u}_{{p}}(t)-\boldsymbol{u}_{{p}}(0){{|}^{2}}\rangle =\frac{1}{{{\kappa}_{p}}}\left(1-{{\text{e}}^{-2t/{{\tau}_{p}}}}\right),~\text{for}~p\geqslant 1, \\ \sigma _{0}^{2}(t) & =2{{D}_{\text{cm}}}t. \end{array}\end{eqnarray} \tag{ 108 }$$

The relaxation times are defined by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\tau}_{p}}=\frac{1}{D{{\kappa}_{p}}}, \end{array}\end{eqnarray} \tag{ 109 }$$

while the diffusion constant is ${{D}_{\text{cm}}}=D/N$. The shortest timescale is ${{\tau}_{N-1}}\approx 1/\left(4D\kappa \right)$ which is half of the time ${{\tau}_{\text{s}}}=1/\left(2D\kappa \right)$ for a free monomer to diffuse a mean squared distance between adjacent monomers (${{b}^{2}}=1/\kappa$). The center of mass is characterized by the time scale ${{\tau}_{0}}\equiv {{b}^{2}}N/{{D}_{\text{cm}}}={{N}^{2}}/\left(D\kappa \right)$ which is the time for a particle to diffusion across the polymer size. For long polymers ${{\tau}_{0}}/{{\tau}_{1}}\approx {{\pi}^{2}}$. Using relation (108), the MSD of monomer $\boldsymbol{R}_{{c}}$ is a sum of independent OU-variables,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{var}\left(\boldsymbol{R}_{{\boldsymbol{c}}}(t)\right)=\langle {{\left(\boldsymbol{R}_{{c}}(t)-\boldsymbol{R}_{{c}}(0)\right)}^{2}}\rangle =\frac{d}{\kappa N}\underset{p=1}{\overset{N-1}{\sum}}\,\frac{{{\cos}^{2}}\left(\frac{(2c-1)p\pi}{2N}\right)}{{{\sin}^{2}}\left(\frac{p\pi}{2N}\right)}\left(1-{{\text{e}}^{-2t/{{\tau}_{p}}}}\right)+2d{{D}_{\text{cm}}}t, \end{array}\end{eqnarray} \tag{ 110 }$$

where d is the spatial dimension. Formula (110) shows the deviation with MSD of a Brownian motion, for which the correlation function increases linearly with time. There are three distinguish regimes:

• 1.  
  For short time $t\ll {{\tau}_{N-1}}$, $\sigma _{p}^{2}(t)\approx Dt$, independent of p, the sum in equation (110) leads to

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{var}\left(\boldsymbol{R}_{{\boldsymbol{c}}}\right)\approx 2dDt, \end{array}\end{eqnarray} \tag{ 111 }$$

  which is the diffusion regime.
• 2.  
  For large time, $t\gg {{\tau}_{1}}$, the exponential terms in relation (110) becomes independent of t. Only the first term in equation (108) corresponding to the diffusion of the center of mass gives the time-dependent behavior. This regime is dominated by normal diffusion, with a diffusion coefficient D/N.
• 3.  
  For intermediate times ${{\tau}_{N-1}}\ll t\ll {{\tau}_{1}}$, such that $2t/{{\tau}_{p}}>1$, the sum of exponentials contributes to equation (110). The variance (110) is

  $$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{var}\left({{R}_{c}}\right)\approx 2{\int}_{{{p}_{\text{min}}}}^{N-1}\frac{{{\cos}^{2}}\left(\frac{(2c-1)p\pi}{2N}\right)}{{{\sin}^{2}}\left(\frac{p\pi}{2N}\right)}\text{d}p, \end{array}\end{eqnarray} \tag{ 112 }$$

  where p_min is such that ${{\tau}_{{{p}_{\min}}}}=2t$. We have $\text{var}\left(\boldsymbol{R}_{{c}}\right)\sim {{t}^{1/2}}$. A Rouse monomer exhibits anomalous diffusion. The time interval can be arbitrarily long with the size N of the polymer.

The first looping time between two monomers is the First Encounter Time ${{\tau}_{e}}$ for two monomers ${{n}_{a}},{{n}_{b}}$ to come into a distance $\varepsilon <b$, defined by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\tau}_{e}}=\inf \left\{t>0~\text{such}~\text{that}~|\boldsymbol{R}_{{{{n}_{a}}}}(t)-\boldsymbol{R}_{{{{n}_{b}}}}(t)|\leqslant \varepsilon \right\}, \end{array}\end{eqnarray} \tag{ 113 }$$

where $\boldsymbol{R}_{{{{n}_{a}}}}$ and $\boldsymbol{R}_{{{{n}_{b}}}}$ follows for example the Rouse equation (98). We present the asymptotic computation of the Mean First Encounter Time (MFET) $\langle {{\tau}_{e}}\rangle$ for the two ends $\boldsymbol{R}_{{N}},\boldsymbol{R}_{{1}}$ meets. The two monomers meet when distance is less than $\varepsilon <b$, that is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} |\boldsymbol{R}_{{N}}-\boldsymbol{R}_{{1}}|\leqslant \varepsilon \end{array}\end{eqnarray} \tag{ 114 }$$

In Rouse coordinates, $\boldsymbol{u}_{{\boldsymbol{p}}}={\sum}_{n=1}^{N}\alpha _{p}^{n}\boldsymbol{R}_{{n}}$ where $\alpha _{p}^{n}$ are defined in (102), condition (114) is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} |2\sqrt{\frac{2}{N}}\underset{p\;\text{odd}}{\sum}\,\boldsymbol{u}_{{p}}\cos \left(\,p\pi /2N\right)|\leqslant \varepsilon. \end{array}\end{eqnarray} \tag{ 115 }$$

The end-to-end encounter is independent of the center of mass, with coordinate $\boldsymbol{u}_{{0}}$. Thus, the MFET is the MFPT for the (N  −  1)d-dimensional stochastic process

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \boldsymbol{u}(t)=\left(\boldsymbol{u}_{{1}}(t),..\boldsymbol{u}_{{N-1}}(t)\right)\in \Omega \times \Omega...\times \Omega =\rm{\tilde \Omega }, \end{array}\end{eqnarray} \tag{ 116 }$$

where $\Omega ={{\mathbb{R}}^{2}}$ or ${{\mathbb{R}}^{3}}$ and $\boldsymbol{u}_{{p}}$ satisfies the OU-equations (105) to the boundary of the domain

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{S}_{\epsilon}}=\left\{P\in \rm{\tilde \Omega }~\text{such}~\text{that}~\text{dist}\left(P,\mathcal{S}\right)\leqslant \frac{\varepsilon}{\sqrt{2}}\right\}, \end{array}\end{eqnarray} \tag{ 117 }$$

where dist is the Euclidean distance and

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \mathcal{S}=\left\{\left(\boldsymbol{u}_{{1}},..\boldsymbol{u}_{{N-1}}\right)\in \rm{\tilde \Omega }|\underset{p\;\text{odd}}{\sum}\,\boldsymbol{u}_{{p}}\cos \left(\,p\pi /2N\right)=0\right\} \end{array}\end{eqnarray} \tag{ 118 }$$

is a submanifold of codimension d in $\rm{\tilde \Omega }$. The probability density function (pdf) $p\left(\boldsymbol{u}(t)=\boldsymbol{x},t\right)$ satisfies the forward Fokker–Planck equation (FPE) [18]

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{1}{D}\frac{\partial p\left(\boldsymbol{x},t\right)}{\partial t} & = & \Delta p\left(\boldsymbol{x},t\right)+\nabla \centerdot \,\left(\nabla \phi \;p\left(\boldsymbol{x},t\right)\right)=\mathcal{L}p, \\ p\left(\boldsymbol{x},0\right) & = & {{p}_{0}}\left(\boldsymbol{x}\right), \end{array}\end{eqnarray} \tag{ 119 }$$

with boundary condition $p\left(\boldsymbol{x},t\right)=0$ for $x\in \partial {{S}_{\epsilon}}$, ${{p}_{0}}\left(\boldsymbol{x}\right)$ is the initial distribution and the potential $\phi \left({{u}_{1}},..,{{u}_{N}}\right)=\frac{1}{2}{\sum}_{p}{{\kappa}_{p}}\boldsymbol{u}_{p}^{2}$ was introduced in (103).

The solution of equation (119) is expanded in eigenfunctions

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} p\left(\boldsymbol{x},t\right)=\underset{i=0}{\overset{\infty}{\sum}}\,{{a}_{i}}{{w}_{\lambda _{i}^{\epsilon}}}\left(\boldsymbol{x}\right){{\text{e}}^{-\lambda _{i}^{\epsilon}tD}}{{\text{e}}^{-\phi \left(\boldsymbol{x}\right)}}, \end{array}\end{eqnarray} \tag{ 120 }$$

where a_i are coefficients, ${{w}_{\lambda _{i}^{\epsilon}}}\left(\boldsymbol{x}\right)$ and $\lambda _{i}^{\epsilon}$ are the eigenfunctions and eigenvalues respectively of the operator $\mathcal{L}$ in the domain ${{ \Omega }_{\epsilon}}=\rm{\tilde \Omega }-{{S}_{\epsilon}}$. The probability distribution that the two ends have not met before time t is the survival probability

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} p(t)=Pr\left\{{{\tau}_{\epsilon}}>t\right\}={{{\int}^{}}_{{{ \Omega }_{\epsilon}}}}p\left(\boldsymbol{x},t\right)\text{d}x, \end{array}\end{eqnarray} \tag{ 121 }$$

and the first looping time is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\tau}_{\epsilon}}=\inf \left\{t>0,\boldsymbol{u}(t)\in \partial {{S}_{\epsilon}}\right\}. \end{array}\end{eqnarray} \tag{ 122 }$$

Using expansion (120), $p(t)=4{\sum}_{i=0}^{\infty}{{C}_{i}}{{\text{e}}^{-\lambda _{i}^{\epsilon}Dt}}$ where ${{C}_{i}}={{{\int}^{}}_{{{ \Omega }_{\epsilon}}}}{{p}_{0}}\left(\boldsymbol{x}\right){{w}_{\lambda _{i}^{\epsilon}}}\left(\boldsymbol{x}\right)\text{d}\boldsymbol{x}{\int}_{{{ \Omega }_{\epsilon}}}{{w}_{\lambda _{i}^{\epsilon}}}\left(\boldsymbol{x}\right){{\text{e}}^{-\phi \left(\boldsymbol{x}\right)}}\text{d}\boldsymbol{x}$. Starting with an equilibrium distribution ${{p}_{0}}\left(\boldsymbol{x}\right)=|\rm{\tilde \Omega }{{|}^{-1}}{{\text{e}}^{-\phi \left(\boldsymbol{x}\right)}}$, we have

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{C}_{i}}=|\rm{\tilde \Omega }{{|}^{-1}}{{\left({{{\int}^{}}_{{{ \Omega }_{\epsilon}}}}{{w}_{\lambda _{i}^{\epsilon}}}\left(\boldsymbol{x}\right){{\text{e}}^{-\phi \left(\boldsymbol{x}\right)}}\text{d}\boldsymbol{x}\right)}^{2}} \end{array}\end{eqnarray}$$

and finally the MFET is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle {{\tau}_{\epsilon}}\rangle =\underset{i=0}{\overset{\infty}{\sum}}\,\frac{{{C}_{i}}}{D\lambda _{i}^{\epsilon}}. \end{array}\end{eqnarray} \tag{ 123 }$$

Starting from the equilibrium distribution ${{C}_{0}}\approx 1$, while the other terms are C_i  =  o(1), as we shall see, the first term is the main contributor of the series.

The eigenvalues $\lambda _{i}^{\epsilon}$ of the operator $\mathcal{L}$ (equation (119)) are obtained by solving the forward FPE in ${{\mathbb{R}}^{d(N-1)}}$, with the zero absorbing boundary condition on the entire boundary of the domain ${{S}_{\epsilon}}$ (see equation (117)), which is the tubular neighborhood of the (N  −  1)d-dimensional sub-manifold $\mathcal{S}$. For small ε, the eigenvalues can be computed from the following regular expansion near the solution of the domain with no domain ${{S}_{\epsilon}}$:

$$\begin{eqnarray}\begin{array}{*{35}{l}} \lambda _{i}^{\epsilon} & = & \lambda _{i}^{0}+{{c}_{2}}\epsilon {\int}_{S}w_{\lambda _{i}^{0}}^{2}\text{d}{{V}_{x}}+\mathcal{O}\left({{\epsilon}^{2}}\right),~\text{for}~d=3 \end{array}\end{eqnarray} \tag{ 124 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \lambda _{i}^{\epsilon} & = & \lambda _{i}^{0}+\frac{2\pi}{\log \epsilon}{{{\int}^{}}_{S}}w_{\lambda _{i}^{0}}^{2}\text{d}{{V}_{x}}+\mathcal{O}\left({{\left(\frac{1}{\log \epsilon}\right)}^{2}}\right)~\text{for}~d=2, \end{array}\end{eqnarray} \tag{ 125 }$$

where the eigenfunction ${{w}_{\lambda _{i}^{0}}}$ and eigenvalues $\lambda _{i}^{0}$ are associated to the non perturbed operator (no boundary) [70], d  =  3,2.

In the context of the Rouse polymer, the volume element is $\text{d}{{V}_{x}}={{\text{e}}^{-\phi \left(\boldsymbol{x}\right)}}\text{d}\boldsymbol{x}_{{g}}$, $\text{d}\boldsymbol{x}_{{g}}$, a measure over the sub-manifold $\mathcal{S}$ and ${{c}_{2}}=\frac{2{{\pi}^{3/2}}}{ \Gamma (3/2)}$ [70]. The unperturbed eigenfunctions ${{w}_{\lambda _{i}^{0}}}$ are products of Hermite polynomials [71], that depend on the spatial coordinates and the eigenvalues $\lambda _{i}^{0}$ are the sum of one dimensional eigenvalues [72]. The first eigenfunction associated to the zero eigenvalue is ${{w}_{\lambda _{0}^{0}}}=|\rm{\tilde \Omega }{{|}^{-1/2}}$. The first eigenvalue for ε small is obtained from relation (124) in dimension 3 with $\lambda _{i}^{0}=0$,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \lambda _{0}^{\epsilon}=\frac{{{c}_{2}}\epsilon {\int}_{S}{{\text{e}}^{-\phi \left(\boldsymbol{x}\right)}}\text{d}\boldsymbol{x}_{{g}}}{|\rm{\tilde \Omega }|}+\mathcal{O}\left({{\epsilon}^{2}}\right){}, \end{array}\end{eqnarray} \tag{ 126 }$$

which is the ratio of the closed to all polymer configurations. Using the potential ϕ (defined in (103)), the volume is computed explicitly from Gaussian integrals

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} |\rm{\tilde \Omega }|={{{\int}^{}}_{ \Omega }}{{\text{e}}^{-\phi \left(\boldsymbol{x}\right)}}\text{d}\boldsymbol{x}_{{g}}={{\left[\frac{{{\left(2\pi \right)}^{(N-1)}}}{\underset{1}{\overset{N-1}{\prod}}\,{{\kappa}_{p}}}\right]}^{d/2}}, \end{array}\end{eqnarray} \tag{ 127 }$$

while the parametrization of the constraint (118) leads to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\int}^{}}_{S}}{{\text{e}}^{-\phi \left(\boldsymbol{x}\right)}}\text{d}\boldsymbol{x}_{{g}}={{\left[\frac{{{\left(2\pi \right)}^{N-2}}\underset{p\;\text{odd}}{\prod}\,\omega _{p}^{2}}{\underset{p}{\prod}\,{{\kappa}_{p}}\left(\underset{p\;\text{odd}}{\sum}\,\frac{\omega _{p}^{2}}{{{\kappa}_{p}}}\right)}\right]}^{d/2}}, \end{array}\end{eqnarray} \tag{ 128 }$$

where ${{\omega}_{p}}=\cos \left(\,p\pi /2N\right)$. In summary, for fix N and small ε,

$$\begin{eqnarray}\lambda _{0}^{\epsilon}=\left\{\begin{array}{*{35}{l}} {{\left(\frac{\kappa}{N\pi}\right)}^{3/2}}4\pi \epsilon +\mathcal{O}\left({{\epsilon}^{2}}\right) & \quad \text{for}~d=3, \\ {} & {} \\ \frac{2\kappa}{N\log \left(\frac{\sqrt{2}b}{\varepsilon}\right)}+\mathcal{O}\left({{\left(\frac{1}{\log \epsilon}\right)}^{2}}\right) & \quad \text{for}~d=2. \end{array}\right. \end{eqnarray} \tag{ 129 }$$

The zero eigenvalue is sufficient to characterize the MFET, confirming that the FET is almost Poissonian, except for very short time. Moreover, the second term in the expansion of $\lambda _{0}^{\epsilon}$ is proportional to 1/N. Using the approximation ${{C}_{0}}\sim 1$ and relation (123), for d  =  3, the MFET is approximated by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\langle {{\tau}_{\varepsilon}}\rangle}_{3d}} & \approx & \frac{1}{D\lambda _{0}^{\epsilon}}=\frac{1}{D\left({{\left(\frac{\kappa}{N\pi}\right)}^{3/2}}4\pi \epsilon -A{{\epsilon}^{2}}/N\right)}, \end{array}\end{eqnarray} \tag{ 130 }$$

where A is a constant that has been estimated numerically. Indeed, with $\epsilon =\frac{\varepsilon}{\sqrt{2}}$, the MFET is for d  =  3

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\langle {{\tau}_{\varepsilon}}\rangle}_{3d}} & = & {{\left(\frac{N\pi}{\kappa}\right)}^{3/2}}\frac{\sqrt{2}}{D4\pi \varepsilon}+{{A}_{3}}\frac{{{b}^{2}}}{D}{{N}^{2}}+\mathcal{O}(1), \end{array}\end{eqnarray} \tag{ 131 }$$

which hold for a large range of N, as evaluated with Brownian simulations (figure 3). The value of the coefficient is A₃  =  0.053 [72] (${{A}_{3}}=0.053{{b}^{2}}/D$). These estimates are obtained for fixed N and small ε.

Zoom In Zoom Out Reset image size

Similar for d  =  2, a a two dimensional space, the asymptotic formula for MFET [72] is

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\langle {{\tau}_{\varepsilon}}\rangle}_{2d}} & = & \frac{N}{2D\kappa}\log \left(\frac{\sqrt{2}b}{\varepsilon}\right)+{{A}_{2}}\frac{{{b}^{2}}}{D}{{N}^{2}}+\mathcal{O}(1), \end{array}\end{eqnarray} \tag{ 132 }$$

All these asymptotic expansions are derived for fix N and small ε. However, there should not be valid in the limit N large, although stochastic simulations (figures 3(a) and (b)) shows that the range validity is broader than expected. The exact asymptotic formula for any two monomers inside a polymer chain should be derived. Other scaling laws have been derived in [68, 69, 73].

A simplified model of a crowded membrane can be a square lattice of circular obstacles of radius a centered at the corners of lattice squares of side L (figure 5). The mean exit time from a lattice box, formula (82), is to leading order independent of the starting position (x,y) and can be approximated as

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\bar{\tau}}}_{4}}=\frac{{\bar{\tau}}}{4}, \end{array}\end{eqnarray} \tag{ 133 }$$

where $\bar{\tau}$ is the MFPT to a single absorbing window in a narrow strait with the other windows closed (reflecting instead of absorbing). It follows that the waiting time in the cell enclosed by the obstacles is exponentially distributed with rate

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \lambda =\frac{1}{\bar{2}{{\tau}_{4}}}, \end{array}\end{eqnarray} \tag{ 134 }$$

where $\bar{\tau}$ is given by (80) and (82) as

$$\begin{eqnarray}\bar{\tau}\approx \left\{\begin{array}{*{35}{l}} {} & {{c}_{1}}~\text{for}~0.8<\varepsilon <1, \\ {} & {} \\ {} & {{c}_{2}}| \Omega |\log \frac{1}{\varepsilon}+{{d}_{1}}~\text{for}~0.55<\varepsilon <0.8, \\ {} & {} \\ {} & {{c}_{3}}\frac{| \Omega |}{\sqrt{\varepsilon}}+{{d}_{2}}~\text{for}~\varepsilon <0.55, \end{array}\right. \end{eqnarray} \tag{ 135 }$$

with $\varepsilon =(L-2a)/a$ and ${{d}_{1}},{{d}_{2}}=O(1)$ for $\varepsilon \ll 1$ (see figure 5). The MFPT c₁ from the center to the boundary of an unrestricted square is computed from

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} u(x,y)=\frac{4{{L}^{2}}}{{{\pi}^{3}}D}\underset{0}{\overset{\infty}{\sum}}\,\frac{\left[\cosh \left(k+\frac{1}{2}\right)\pi -\cosh \left(k+\frac{1}{2}\right)\pi \left(2y/L-1\right)\right]\sin (2k+1)\pi x/L}{{{(2k+1)}^{3}}\cosh (2k+1)\pi}, \end{array}\end{eqnarray} \tag{ 136 }$$

so ${{c}_{1}}=u\left(L/2,L/2\right)\approx \left[4{{L}^{2}}/{{\pi}^{3}}D\right]\left[\cosh \left(\pi /2-1\right)/\cosh \pi \right].$ For L  =  1, D  =  1, we find ${{c}_{1}}\approx 0.076\,$, in agreement with Brownian dynamics simulations (figure 5(B)). The coefficient c₂ is obtained from (78) as ${{c}_{2}}=1/2\pi D\approx 0.16.$ Similarly, the coefficient c₃ is obtained from (82) as ${{c}_{3}}\approx \pi /4\sqrt{2}\,D\approx 0.56.$ The coefficients d_i are chosen by patching $\bar{\tau}$ continuously between the different regimes:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{d}_{1}}={{c}_{1}}+{{c}_{2}}| \Omega \left({{r}_{1}}\right)|\log \left(1-2{{r}_{1}}\right), \end{array}\end{eqnarray} \tag{ 137 }$$

and

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{d}_{2}}={{c}_{1}}-{{c}_{2}}\left[| \Omega \left({{r}_{1}}\right)|\log \left(1-2{{r}_{1}}\right)\,+\,| \Omega \left({{r}_{2}}\right)|\log \left(1-2{{r}_{2}}\right)\right] \\ \,\,\,\,\,\,\,\,\,\,\,\,-\,\,{{c}_{3}}| \Omega \left({{r}_{2}}\right)|{{\left(1-2{{r}_{2}}\right)}^{-1/2}}, \end{array}\end{eqnarray}$$

where $| \Omega (r)|={{L}^{2}}-\pi {{r}^{2}}$.

Zoom In Zoom Out Reset image size

Simulations with D  =  1 in a square of radius L  =  1 with four reflecting circles of radius r, centered at the corners, show that the uniform approximation by the patched formula (135) is in good agreement with Brownian results (figure 5(b)), where the statistics were collected from 1000 escape times of Brownian trajectories per graph point. The trajectories start at the square center. Equation (135) holds in the full range of values of $a\in \left[0,L/2\right]$ and all L.

The Brownian motion around the obstacles (figure 5(a)) can be coarse-grained into a Markovian jump process whose state are the connected domains enclosed by the obstacles and the jump rates are determined from the reciprocals of the mean first passage times and exit probabilities. This random walk can in turn be approximated by an effective coarse-grained anisotropic diffusion. The diffusion approximation to the transition probability density function of an isotropic random walk that jumps at exponentially distributed waiting times with rate λ on a square lattice with step size L is given by [18]

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\partial p}{\partial t}=\bar{D}\left(\frac{{{\partial}^{2}}p}{\partial {{x}^{2}}}+\frac{{{\partial}^{2}}p}{\partial {{y}^{2}}}\right),\quad \bar{D}=\frac{\lambda {{L}^{2}}}{4}. \end{array}\end{eqnarray} \tag{ 138 }$$

The results of the previous section can be used to estimate the density of obstacles on the membrane of cell such as a neuronal dendrite. The effective diffusion coefficient of a receptor on the neuronal membrane can be estimated from the experimentally measured single receptor trajectory by a single particle tracking method. The receptor effective diffusion coefficient of a receptor varies between 0.01 and 0.2 μm²/sec.

In the simplified model of crowding, the circular obstacles are as in (figure 5). Simulated Brownian trajectories give the MFPT from one square to the next one as shown in figure 4, where L is fixed and a is variable. According to (135), (134) and (138), as a increases the effective diffusion coefficient $\bar{D}$ decreases. It is computed as the as the mean square displacement (MSD) $\langle \frac{MSD(t)}{4t}\rangle$. Brownian simulations show that $\bar{D}$ is linear, thus confirming that in the given geometry crowding does not affect the nature of the Brownian motion for sufficiently long times. Specifically, for Brownian diffusion coefficient D  =  0.2 μm² s⁻¹ the time considered is longer than 10 s. In addition, figure 6(c) shows the diffusion coefficient ratio ${{D}_{a}}/{{D}_{0}}$, where D_a is the effective diffusion coefficient of Brownian motion on the square lattice described above with obstacles of radius a. For a  =  0.3 the value ${{D}_{a}}/{{D}_{0}}\approx 0.7$ is found whereas a direct computation from the mean exit time formula (135) gives

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{\tau}_{0}}}{{{\tau}_{a}}}=\frac{{{c}_{1}}}{{{c}_{2}}| \Omega |\log \frac{1}{\varepsilon}+{{d}_{1}}}\approx 0.69, \end{array}\end{eqnarray} \tag{ 139 }$$

where $\varepsilon =(L-2a)/L=0.4$.

Zoom In Zoom Out Reset image size

It can be concluded from the Brownian simulations that the coarse-grained motion is plain diffusion with effective diffusion coefficient ${{D}_{a}}/{{D}_{0}}={{\tau}_{0}}/{{\tau}_{a}}$, which decreases nonlinearly as a function of the radius a, as given by the uniform formula (135). Figure 6 recovers the three regimes of (135): the uncrowded regime for a  <  0.2L, where the effective diffusion coefficient does not show any significant decrease, a region 0.2L  <  a  <  0.4L, where the leading order term of the effective diffusion coefficient is logarithmic, and for a  >  0.4L the effective diffusion coefficient decays as $\sqrt{(L-2a)/L}$, in agreement with (135).

Finally, to estimate the density of obstacles in a neuron from (135), (134) and (138), a reference density has to be chosen. The reference diffusion coefficient is chosen to be that of receptors moving on a free membrane (with removed cholesterol), estimated to be $0.17\leqslant D\leqslant 0.2$ μm² s⁻¹ [14], while with removing actin, the diffusion coefficient is 0.19 μm² s⁻¹. The reference value D  =  .2 μm² s⁻¹ gives an estimate of the crowding effect based on the measured diffusion coefficient (figure 6(d)). The reduction of the diffusion coefficient from D  =  0.2 μm² s⁻¹ to D  =  0.04 μm² s⁻¹ is achieved when 70% of the membrane surface is occupied by obstacles. Thus obstacles impair the diffusion of receptors and are therefore responsible for the large decrease of the measured diffusion coefficient (up to five times). To conclude, as illustrated in figure 7, diffusion in a crowded membrane involves various time regimes: at very short time scale, before a Brownian particle has the time to escape between small aperture near obstacles, the particle diffuses freely, characterized by the homogeneous membrane diffusion coefficient. This approximation is valid before the NET time scale is reached (see formula (82)). For longer time, $t\gg \bar{\tau}$, the motion of a Brownian particle is characterized again by a diffusion process, but now the diffusion coefficient does account for the obstacles and in the very density limit, the effective diffusion coefficient is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{D}_{e}}=D\frac{{{L}^{2}}}{2\pi \left({{L}^{2}}-\pi {{R}^{2}}\right)}\sqrt{\frac{L-2R}{R}},~\text{for}~R\leqslant L/2. \end{array}\end{eqnarray} \tag{ 140 }$$

Zoom In Zoom Out Reset image size

At an intermediate time regime between the two extreme cases described above, a stochastic particle is hopping from one square to another, characterized as anomalous diffusion (blue curves in figure 7).

The model of the telomere dynamics is

$$\begin{eqnarray}{{x}_{n+1}}=\left\{\begin{array}{*{35}{l}} {{x}_{n}}-a & \text{w}.\text{p}.~l\left({{x}_{n}}\right) \\ {} & {} \\ {{x}_{n}}+b & \text{w}.\text{p}.~r\left({{x}_{n}}\right), \end{array}\right. \end{eqnarray} \tag{ 141 }$$

where the right-probability r(x) can be approximated by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} r(x)=\frac{1}{1+\alpha x}, \end{array}\end{eqnarray} \tag{ 142 }$$

for some $\alpha >0$. Scaling ${{x}_{n}}={{y}_{n}}L$ and setting $\varepsilon =b/L$, the dynamics (141) becomes

$$\begin{eqnarray}{{y}_{n+1}}=\left\{\begin{array}{*{35}{l}} {{y}_{n}}-\frac{\varepsilon a}{b} & \text{w}.\text{p}.~\tilde{l}\left(\,{{y}_{n}}\right) \\ {{y}_{n}}+\varepsilon & \text{w}.\text{p}.~\tilde{r}\left(\,{{y}_{n}}\right), \end{array}\right. \end{eqnarray} \tag{ 143 }$$

where $\tilde{l}\left(\,y\right)=l(x)$ and $\varepsilon T/b<y<1$. In the limit $\varepsilon \ll 1$, the process y_n moves in small steps. The dynamics (143) falls under the general scheme (see [18])

$$\begin{eqnarray}&&{{y}_{n+1}}={{y}_{n}}+\varepsilon {{\xi}_{n}},\end{eqnarray} \tag{ 144 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \text{Pr}\left\{{{\xi}_{n}}=\xi \,|\,{{y}_{n}}=y,\,{{y}_{n-1}}={{y}_{1}},\,\ldots \right\}=w\left(\xi \,|\,y,\varepsilon \right), \end{array}\end{eqnarray} \tag{ 145 }$$

ε is a small parameter, and y₀ is a random variable with a given pdf p₀(y). In the case at hand the function $w\left(\xi \,|\,y\right)$ defined in (145) is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} w\left(\xi \,|\,y\right)=\left(1-\tilde{{r}}\,\left(\,y\right)\right)\delta \left(\xi +\frac{a}{b}\right)+\tilde{{r}}\,\left(\,y\right)\delta \left(\xi -1\right). \end{array}\end{eqnarray} \tag{ 146 }$$

The pdf of y_n satisfies the backward equation for 0  <  x  <  1

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & \,{{p}_{\varepsilon}}\left(\,y,n\,|\,x,m\right)-{{p}_{\varepsilon}}\left(\,y,n\,|\,x,m+1\right) \\ = & \,{{p}_{\varepsilon}}\left(\,y,n\,|\,x-\varepsilon \frac{a}{b},m+1\right)\tilde{{l}}\,(x)+{{p}_{\varepsilon}}\left(\,y,n\,|\,x+\varepsilon,m+1\right)\tilde{{r}}\,(x)-{{p}_{\varepsilon}}\left(\,y,n\,|\,x,m\right). \end{array}\end{eqnarray} \tag{ 147 }$$

The first conditional jump moment, ${{m}_{1}}\left(\,y\right)=-a\tilde{l}\left(\,y\right)/b+\tilde{r}\left(\,y\right)$ changes sign at ${{z}_{0}}=b/L\alpha a$, so ${{p}_{\varepsilon}}\left(\,y,n\right)$ converges to a quasi-stationary density ${{p}_{\varepsilon}}\left(\,y\right)$ for large n, before the trajectory y_n is terminated at y  =  T/L.

One dimensional processes (see equation (144)) in the small jump limit are described in [18], p 236 and p 303, see also [76–78]. In this limit, the Kramers–Moyal approximation consists in expanding in ε and then approximating equation (147) by a direct truncation to a second order equation. The method is equivalent to construct a stochastic processes with the first and second moments that match the coefficients of the Kramers–Moyal approximation. We present below the WKB construction of an approximated solution.

The structure of the quasi-stationary solution ${{p}_{\varepsilon}}\left(\,y\right)$ for $\varepsilon \ll 1$ and $y>\varepsilon$, can be obtained in the WKB form [18]

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{\varepsilon}}\left(\,y\right)={{K}_{\varepsilon}}\left(\,y\right)\exp \left\{-\frac{\psi \left(\,y\right)}{\varepsilon}\right\}, \end{array}\end{eqnarray} \tag{ 148 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{K}_{\varepsilon}}\left(\,y\right)=\underset{i=0}{\overset{\infty}{\sum}}\,{{\varepsilon}^{i}}{{K}_{i}}\left(\,y\right) \end{array}\end{eqnarray} \tag{ 149 }$$

and K_i(y) are chosen such that ${{p}_{\varepsilon}}\left(\,y\right)$ is normalized. The components of the expansion (148) are found by solving the eikonal equation for $\psi \left(\,y\right)$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\int}^{}}_{\mathbb{R}}}\left[{{\text{e}}^{\xi {{\psi}^{\prime}}\left(\,y\right)}}-1\right]w\left(\xi \,|\,y,0\right)\,\text{d}\xi =0, \end{array}\end{eqnarray} \tag{ 150 }$$

which for the problem at hand takes the form

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tilde{l}\left(\,y\right)\exp \left\{-{{\psi}^{\prime}}\left(\,y\right)\frac{a}{b}\right\}+\tilde{r}\left(\,y\right)\exp \left\{{{\psi}^{\prime}}\left(\,y\right)\right\}=1~\text{for}~y>\varepsilon. \end{array}\end{eqnarray} \tag{ 151 }$$

At the next order in ε, we find that K₀(y) is the solution of the 'transport' equation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\int}^{}}_{\mathbb{R}}}\left\{\frac{\partial}{\partial y}\left[w\left(\xi \,|\,y\right){{K}_{0}}\left(\,y\right)\right]+\frac{\xi w\left(\xi \,|\,y\right)}{2}{{\psi}^{\prime \prime}}\left(\,y\right){{K}_{0}}\left(\,y\right)\right\}\xi {{\text{e}}^{\xi {{\psi}^{\prime}}\left(\,y\right)}}\,\text{d}\xi =0. \end{array}\end{eqnarray} \tag{ 152 }$$

Equation (152) has a removable singularity at y  =  z₀, because the coefficient of K_0y in (152),

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\int}^{}}_{\mathbb{R}}}\xi w\left(\xi \,|\,y\right){{\text{e}}^{\xi {{\psi}^{\prime}}\left(\,y\right)}}\,\text{d}\xi & = \,{\int}_{\mathbb{R}}\xi w\left(\xi \,|\,y\right)\left[1+\xi {{\psi}^{\prime \prime}}\left({{z}_{0}}\right)\left(\,y-{{z}_{0}}\right)+O\left({{\left(\,y-{{z}_{0}}\right)}^{2}}\right)\right]\,\text{d}\xi \\ & = \left[m_{1}^{\prime}\left({{z}_{0}}\right)+{{m}_{2}}\left({{z}_{0}}\right){{\psi}^{\prime \prime}}\left({{z}_{0}}\right)\right]\left(\,y-{{z}_{0}}\right)+O\left({{\left(\,y-{{z}_{0}}\right)}^{2}}\right) \\ & = -m_{1}^{\prime}\left({{z}_{0}}\right)\left(\,y-{{z}_{0}}\right)+O\left({{\left(\,y-{{z}_{0}}\right)}^{2}}\right), \end{array}\end{eqnarray}$$

vanishes linearly at z₀. However, the coefficient of K₀(y) in (152) also vanishes at z₀.

In this section, we present a different approach for computing the steady solution of the pdf for the process (141). We derive the Takacs equation and then study the statistics of the shortest trajectory (shortest telomere) for an ensemble of n identical independently distributed (iid) processes.

The steady state distribution associated to the telomere equation (141) can be rescaled with a constant drift for shortening, and possible large jumps with exponential rates for elongation. The jump rate function becomes $\lambda (X)=\frac{1}{1+BX}$, and the probability for the jump $\bar{\xi}$ is given by $Pr\left(\bar{\xi}=y\right)=\theta {{\text{e}}^{-\theta y}}=b\left(\,y\right)$, where $B=a\beta$ and $\theta =ap$. The pdf $f(x,t)=\partial F(x,t)/\partial x$ where $F(x,t)=\text{Pr}\left\{X(t)\leqslant x\right\}$ satisfies the Takacs equation, which is written for the forward Fokker–Planck equation x  >  0,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\partial}_{t}}f={{\partial}_{x}}\,f-\lambda (x)f(x,t)+{\int}_{0}^{l}\lambda \left(\,y\right)f\left(\,y,t\right)b(x-y)\text{d}y. \end{array}\end{eqnarray} \tag{ 153 }$$

The stationary distribution function is [75]

$$\begin{eqnarray}\begin{array}{*{35}{l}} \bar{f}(x) & = & \frac{\theta {{\left[\theta \left(x+\frac{1}{B}\right)\right]}^{\frac{1}{B}}}{{\text{e}}^{-\theta \left(x+\frac{1}{B}\right)}}}{ \Gamma \left(\frac{1}{B}+1,\frac{\theta}{B}\right)}, \end{array}\end{eqnarray} \tag{ 154 }$$

where $\Gamma (s,x)={\int}_{x}^{+\infty}{{t}^{s-1}}{{\text{e}}^{-t}}\text{d}t$ is the upper incomplete Gamma function. Now the distribution of the shortest telomere [75] in an ensemble of 2n telomeres, corresponding to a total of n chromosomes (16 in yeast and n is in the range of 36–60) is estimated when their lengths are independent identically distributed variables ${{L}_{1}},{{L}_{2}},\ldots {{L}_{2n}}$. Considering 2n iid variables ${{X}_{1}},\ldots {{X}_{2n}}$ following a distribution f, the pdf of the minimum ${{X}_{(1:2n)}}=\min \left({{X}_{1}},{{X}_{2}},\ldots {{X}_{2n}}\right)$ is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{f}_{{{X}_{(1:2n)}}}}(x)=2n{{(1-F)}^{2n-1}}(x)f(x), \end{array}\end{eqnarray} \tag{ 155 }$$

where $F(x)={\int}_{0}^{x}f(u)du$. The statistical moments $\bar{X}_{(1:2n)}^{k}={{{\int}^{}}_{\mathbb{R}\text{+}}}{{x}^{k}}{{f}_{{{X}_{(1:2n)}}}}(x)\text{d}x$ are given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{X}_{(1:2n)}^{k}=k{{{\int}^{}}_{\mathbb{R}\text{+}}}{{x}^{k-1}}{{(1-F)}^{2n}}(x)\text{d}x. \end{array}\end{eqnarray} \tag{ 156 }$$

When f is a Gamma distribution of parameter α and n sufficiently large, $\bar{X}_{(1:2n)}^{k}$ can be estimated using the Laplace's method.

In the limit x tends to 0, equation (154) with $\alpha =1+\frac{1}{B}$ satisfies $F(x)\approx m{{x}^{r}}$ with $m=\frac{1}{\alpha \Gamma \left(\alpha \right)}>0$ and $r=\alpha >1$,

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{X}_{(1:2n)}^{k}=\frac{k}{r}{\int}_{0}^{+\infty}{{x}^{k/r-1}}\exp \left[2n\ln \left(1-F\left({{x}^{1/r}}\right)\right)\right]\text{d}x\approx \bar{X}_{(1:2n)}^{k}\approx \frac{k \Gamma \left(\frac{k}{\alpha}\right)}{\alpha}{{\left(\frac{\alpha \Gamma \left(\alpha \right)}{2n}\right)}^{k/\alpha}}. \end{array}\end{eqnarray} \tag{ 157 }$$

Using formulae (155)–(157), the pdf and the moments of the shortest telomere length L_1:2n for k  =  1 can be estimated and the shortest telomere length is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\bar{L}}}_{(1:2n)}}\approx {{L}_{0}}-\frac{1}{B}+\frac{ \Gamma \left(1+\frac{1}{B}\right) \Gamma {{\left(\frac{1}{1+1/B}\right)}^{\frac{1}{1+1/B}}}}{p{{\left(1+\frac{1}{B}\right)}^{\frac{1}{1+B}}}{{(2n)}^{\frac{1}{1+1/B}}}}. \end{array}\end{eqnarray} \tag{ 158 }$$

Using the values p  =  0.026, $\beta =0.045$ (B  =  0.16  <  0.5) and L₀  =  90 (yeast) and equation (156) for k  =  1 and 2, the mean shortest telomere length is $184\pm 25$ bps.

To estimate the gap between the shortest telomere and the others, we shall compute the distribution of the second shortest length X_(2:2n). The pdf ${{f}_{{{X}_{(2:2n)}}}}$ of X_(2:2n) is given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{f}_{{{X}_{(2:2n)}}}}(x)=2n(2n-1)F(x){{\left(1-F(x)\right)}^{2n-2}}f(x), \end{array}\end{eqnarray} \tag{ 159 }$$

and the statistical moments $\bar{X}_{(2:2n)}^{k}$ satisfy the induction relation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \bar{X}_{(2:2n)}^{k}=n\bar{X}_{\left(1:2n-1\right)}^{k}-(2n-1)\bar{X}_{(1:2n)}^{k}. \end{array}\end{eqnarray} \tag{ 160 }$$

Using equation (157) for k  =  1, we obtain that the ratio $\frac{{{{\bar{X}}}_{(2:2n)}}}{{{{\bar{X}}}_{(1:2n)}}}$, for n or B  ≫  1 is given asymptotically for $\alpha =1+\frac{1}{B}$ by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{{{{\bar{X}}}_{(2:2n)}}}{{{{\bar{X}}}_{(1:2n)}}}\approx 1+\frac{1}{\alpha}. \end{array}\end{eqnarray} \tag{ 161 }$$

To conclude, for a pdf with a nonzero first order derivative at 0, this ratio is a universal number $\frac{3}{2}$. In the case of yeast, equation (160) reveals that the mean length of the second shortest telomere is 207 bps. Thus, the shortest telomere is on average 22 bps shorter than the second one. This gap results from the statistical property of the telomere number and dynamics and should exist in all species. It suggests that the length of the shortest telomere controls the number of division. The computation of the mean time to threshold is more involved as it requires finding an approximation of the pdf in two time intervals. Interestingly, this time depend on the escape of a coarse-grained stochastic dynamics from an effective potential [81].

To compute the number of proteins, we need to follow simultaneously three variables: the state s of element A, the number m of Krox20 mRNA and the number n of unbound Krox20 proteins. The joint probability p_s(m,n,t) to find element A in state s with m Krox20 mRNA molecules and n free Krox20 proteins is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{s}}\left(m,n,t\right)=Pr\left\{s(t)=s,mRNA(t)=m,Krox20(t)=n\right\} \end{array}\end{eqnarray} \tag{ 162 }$$

and it satisfies a Master equation [18], where only one binding or unbinding occurs at a time

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\partial {{p}_{s}}}{\partial t}\left(m,n,t\right)={{ \Phi }_{s}}{{p}_{s}}\left(m-1,n,t\right)+(m+1) \Psi {{p}_{s}}\left(m+1,n,t\right)-\left({{ \Phi }_{s}}+m \Psi \right){{p}_{s}}\left(m,n,t\right) \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\,\,m\phi {{p}_{s}}\left(m,n-1,t\right)+(n+1)\psi {{p}_{s}}\left(m,n+1,t\right)-\left(m\phi +n\psi \right){{p}_{s}}\left(m,n,t\right) \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\,\,{{\mu}_{s+1}}{{p}_{s+1}}\left(m,n-1,t\right)+(n+1){{\lambda}_{s-1}}{{p}_{s-1}}\left(m,n+1,t\right)-\left({{\mu}_{s}}+n{{\lambda}_{s}}\right){{p}_{s}}\left(m,n,t\right). \end{array}\end{eqnarray} \tag{ 163 }$$

The mRNA production rate is given function

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{ \Phi }_{s}}(t)={{ \Phi }_{I}}\theta \left({{t}_{I}}-t\right)+{{ \Phi }_{A,s}} \end{array}\end{eqnarray} \tag{ 164 }$$

where $\theta (t)$ is the Heaviside function. There function models an initial phase where proteins are produced until time t_I. During this phase, mRNA is produced with a Poissonian rate ${{ \Phi }_{I}}$ by an external molecule. Each mRNA protein is degraded with a Poissonian rate $\Psi$. Proteins are produced with a Poissonian rate ϕ and are degraded with a rate ψ and can bind to a promoter site that has N_b  =  4 binding sites. The state of A is characterized by s  =  0, 1, 2, 3, 4 bound molecules. The autoregulatory production of mRNA occurs with Poissonian rates ${{ \Phi }_{A,s}}={{ \Phi }_{A}}{{\xi}_{s}}$ that depend on state s, where ${{ \Phi }_{A}}$ is the maximal production rate and ${{\xi}_{s}}$ describes the modulation due to the state of element A. Binding and unbinding of proteins to A are described by state dependent binding and unbinding rates ${{\lambda}_{s}}$ and ${{\mu}_{s}}$.

When the change in the number of free Krox20 proteins n due to binding and unbinding to element A is neglected, in the limit $n\gg {{N}_{b}}$, the Master equation (163) can be approximated by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \frac{\partial {{p}_{s}}}{\partial t}\left(m,n,t\right)={{ \Phi }_{s}}{{p}_{s}}\left(m-1,n,t\right)+(m+1) \Psi {{p}_{s}}\left(m+1,n,t\right)-\left({{ \Phi }_{s}}+m \Psi \right){{p}_{s}}\left(m,n,t\right) \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\,\,m\phi {{p}_{s}}\left(m,n-1,t\right)+(n+1)\psi {{p}_{s}}\left(m,n+1,t\right)-\left(m\phi +n\psi \right){{p}_{s}}\left(m,n,t\right) \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+\,\,{{\mu}_{s+1}}{{p}_{s+1}}\left(m,n,t\right)+n{{\lambda}_{s-1}}{{p}_{s-1}}\left(m,n,t\right)-\left({{\mu}_{s}}+n{{\lambda}_{s}}\right){{p}_{s}}\left(m,n,t\right). \end{array}\end{eqnarray} \tag{ 165 }$$

The first line in equation (165) describes the production and degradation of a mRNA, while the second is for the production and degradation of a Krox20 protein. The last one is for the binding and unbinding of a Krox20 protein to element A. The marginal probabilities is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{s}}\left(m,n,t\right)=\underset{{{s}^{\prime}}=0}{\overset{{{N}_{b}}}{\sum}}\,{{p}_{s,{{s}^{\prime}}}}\left(m,n,t\right)=\underset{{{s}^{\prime}}=0}{\overset{{{N}_{b}}}{\sum}}\,{{p}_{{{s}^{\prime}},s}}\left(m,n,t\right). \end{array}\end{eqnarray} \tag{ 166 }$$

Binding and unbinding to element A is fast compared to the turn over of proteins, thus we use the approximation

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{s}}\left(m,n,t\right)\approx \,p\left(m,n,t\right){{p}_{s}}(n), \end{array}\end{eqnarray} \tag{ 167 }$$

where p_s(n) are the steady-state probabilities to find element A in state s for a given number of Krox20 proteins n, and p(m,n,t) is the probability to find m mRNA molecules and n proteins at time t.

The steady state condition for binding and unbinding from equation (165) is

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mu}_{s+1\,}}{{p}_{s+1}}(n)=n{{\lambda}_{s}}{{p}_{s}}(n), \end{array}\end{eqnarray} \tag{ 168 }$$

leading to the solution

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{p}_{s}}(n)=\frac{\underset{j=0}{\overset{s-1}{\prod}}\,n{{\gamma}_{j}}}{\underset{i=0}{\overset{{{N}_{b}}}{\sum}}\,\underset{j=0}{\overset{i-1}{\prod}}\,n{{\gamma}_{j}}}, \end{array}\end{eqnarray} \tag{ 169 }$$

where

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\gamma}_{i}}=\frac{{{\lambda}_{i}}}{{{\mu}_{i+1}}},\quad i=0,\ldots {{N}_{b}}-1. \end{array}\end{eqnarray} \tag{ 170 }$$

The effective mRNA production rates is defined by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{ \Phi }_{A}}(n)=\underset{s=0}{\overset{{{N}_{b}}}{\sum}}\,{{ \Phi }_{A,s}}{{p}_{s}}(n)={{ \Phi }_{A}}\underset{s=0}{\overset{{{N}_{b}}}{\sum}}\,{{\xi}_{s}}{{p}_{s}}(n). \end{array}\end{eqnarray} \tag{ 171 }$$

At this stage, the approximated Master equation for the joint pdf p(m,n,t) is

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{\partial {{p}_{s}}}{\partial t}\left(m,n,t\right)= & = & \Phi (n,t)p\left(m-1,n,t\right)+(m+1) \Psi p\left(m+1,n,t\right)-\left(\Phi (n,t)+m \Psi \right)p\left(m,n,t\right) \\ {} & {} & +\,\,m\phi p\left(m,n-1,t\right)+(n+1)\psi p\left(m,n+1,t\right)-\left(m\phi +n\psi \right)p\left(m,n,t\right). \end{array}\end{eqnarray} \tag{ 172 }$$

The marginal probabilities for mRNA molecules, Krox20 proteins and the state of element A are

$$\begin{eqnarray}\begin{array}{*{35}{l}} p(m,t) & = & \underset{s=0}{\overset{{{N}_{b}}}{\sum}}\,\underset{n=0}{\overset{\infty}{\sum}}\,{{p}_{s}}\left(m,n,t\right)=\underset{n=0}{\overset{\infty}{\sum}}\,p\left(m,n,t\right), \\ p(n,t) & = & \underset{s=0}{\overset{{{N}_{b}}}{\sum}}\,\underset{m=0}{\overset{\infty}{\sum}}\,{{p}_{s}}\left(m,n,t\right)=\underset{m=0}{\overset{\infty}{\sum}}\,p\left(m,n,t\right), \\ {{p}_{s}}(t) & = & \underset{m=0}{\overset{\infty}{\sum}}\,\underset{n=0}{\overset{\infty}{\sum}}\,{{p}_{s}}\left(m,n,t\right)=\underset{n=0}{\overset{\infty}{\sum}}\,p(n,t){{p}_{s}}(n). \end{array}\end{eqnarray} \tag{ 173 }$$

It remains difficult to study equations (172) because n can be large while m is of the order of few and thus there is no clear limit approximations. The long-time asymptotic is however an important quantity to estimate and in particular how it depends on initial conditions. This limit tells us whether or not a cell expresses Krox20. The diffusion approximation leads to a dynamical system that predicts a two state attractors characterized by full or zero expression, while the Master system shows that there is be a continuum level of expression, but the distribution is dominated by two peaks. This difference justifies the need of a Markov chain description, in particular to study cells at the boundary between regions, where a graded expression is predicted.

The mean field approximation for the mean quantities is based on the scaled variables

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \tau = \Psi t,\quad \hat{m}=\frac{m}{{{m}_{0}}},\quad \hat{n}=\frac{n}{{{n}_{0}}}~ \end{array}\end{eqnarray} \tag{ 174 }$$

and the normalized parameters are

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{m}_{0}}=\frac{{{ \Phi }_{A}}}{ \Psi },\quad {{n}_{0}}={{m}_{0}}\frac{\phi}{\psi},\quad \alpha =\beta {{n}_{0}},\quad \epsilon =\frac{ \Psi }{\psi},\quad \chi =\frac{{{ \Phi }_{I}}}{{{ \Phi }_{A}}}, \end{array}\end{eqnarray} \tag{ 175 }$$

where the average values m₀ and n₀ characterize the mRNA and proteins when element A is fully activated. The mean-field equation in the scaled variables $\hat{m}$ and $\hat{n}$ is expressed using the function

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} f(x)=\underset{s=0}{\overset{{{N}_{b}}}{\sum}}\,{{\xi}_{s}}\frac{\underset{j=0}{\overset{s-1}{\prod}}\,x{{{\tilde{\gamma}}}_{j}}}{\underset{i=0}{\overset{{{N}_{b}}}{\sum}}\,\underset{j=0}{\overset{i-1}{\prod}}\,x{{{\tilde{\gamma}}}_{j}}}. \end{array}\end{eqnarray} \tag{ 176 }$$

Indeed, the Kramers–Moyal expansion of equation (172) is

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{\partial p\left(\hat{m},\hat{n},\tau \right)}{\partial \tau} & = & \underset{k=1}{\overset{\infty}{\sum}}\,{{\left(\frac{1}{{{m}_{0}}}\right)}^{k-1}}\frac{\partial _{{\hat{m}}}^{k}}{k!}\left({{\left(-1\right)}^{k}}\left(\chi \theta \left({{\tau}_{I}}-\tau \right)+f\left(\alpha \hat{n}\right)\right)+\hat{m}\right)p\left(\hat{m},\hat{n},\tau \right) \\ {} & {} & +\,\,\frac{1}{\epsilon}\underset{k=1}{\overset{\infty}{\sum}}\,{{\left(\frac{1}{{{n}_{0}}}\right)}^{k-1}}\frac{\partial _{{\hat{n}}}^{k}}{k!}\left({{\left(-1\right)}^{k}}\hat{m}+\hat{n}\right)p\left(\hat{m},\hat{n},\tau \right). \end{array}\end{eqnarray} \tag{ 177 }$$

From truncating the series at first order and neglecting the initiation ($\chi =0$), we obtain the first order dynamical system, which is the mean-field equation [18]:

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{\text{d}}{\text{d}\tau}\hat{m} & = & \chi \theta \left({{\tau}_{I}}-\tau \right)+f\left(\alpha \hat{n}\right)-\hat{m} \\ {} & {} & {} \\ \frac{\text{d}}{\text{d}\tau}\hat{n} & = & \frac{1}{\epsilon}\left(\hat{m}-\hat{n}\right). \end{array}\end{eqnarray} \tag{ 178 }$$

The fixed points are given by ${{\hat{n}}^{\ast}}={{\hat{m}}^{\ast}}$ in equation (178). When $\chi =0$, there is a minimal value ${{\alpha}_{\min}}\approx 33.8$ for which $\alpha <{{\alpha}_{\min}}$, there is a single fixed point ${{\hat{m}}^{\ast}}=0$ (${{\alpha}_{\min}}$, computed using the solution z₀ of $\frac{{{f}^{\prime}}(x)}{f(x)}=\frac{1}{x}$. The minimum value is ${{\alpha}_{\min}}=\frac{{{z}_{0}}}{f\left({{z}_{0}}\right)}=\frac{1}{{{f}^{\prime}}\left({{z}_{0}}\right)}$).

For $\alpha >{{\alpha}_{\min}}$, there are two stable fixed points $\hat{m}_{1}^{\ast}$ and $\hat{m}_{3}^{\ast}$ and an unstable saddle point $\hat{m}_{2}^{\ast}$ (figure 9 upper). For large α, the asymptotic values are $\hat{m}_{2}^{\ast}\to 0$ and $\hat{m}_{3}^{\ast}\to 1$. The value $\hat{m}_{3}^{\ast}\approx 1$ ($m_{3}^{\ast}={{m}_{0}}\hat{m}_{3}^{\ast}\approx \frac{{{\phi}_{A}}}{ \Psi }$) corresponds to the situation where element A is fully activated (figure 9).

Zoom In Zoom Out Reset image size

The two stable fixed points defines two basins of attraction and a saddle point for $\alpha >{{\alpha}_{\min}}$, contained in a separatrix (figure 9 upper for $\alpha =40$). For an initial conditions outside the basin of attraction of the fixed point $\hat{m}_{1}^{\ast}$, the dynamics evolves towards the Up attractor (a high protein expression level). In contrast, for a small amplitude initiation, the expression vanishes. In figure 9 (upper left) shows the separatrix for different values of α: with increasing α the basin of attraction of $\hat{m}_{1}^{\ast}=0$ shrinks and asymptotically vanishes for large α. In summary, the mean field dynamical system underlying protein activation shows a bistable behavior between two attractors depending on the initial condition, which can be seen as a random variable In contrast, the numerical analysis of the Master equation (172) for the probability p(m,n,t) with the scaled variables $\hat{m}$ and $\hat{n}$, $p\left(\hat{m},\hat{n},t\right)$ for p(m,n,t) uses the marginal probabilities defined in equation (173). The mean values for $\hat{m}$ and $\hat{n}$ are given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} \overline{{\hat{m}}}(t) & = & \frac{1}{{{m}_{0}}}\underset{n=0}{\overset{\infty}{\sum}}\,mp\left(m,n,t\right)=\frac{1}{{{m}_{0}}}\underset{m=0}{\overset{\infty}{\sum}}\,mp(m,t)=\frac{\bar{m}(t)}{{{m}_{0}}}, \\ \overline{{\hat{n}}}(t) & = & \frac{1}{{{n}_{0}}}\underset{m,n=0}{\overset{\infty}{\sum}}\,np\left(m,n,t\right)=\frac{1}{{{n}_{0}}}\underset{n=0}{\overset{\infty}{\sum}}\,np(n,t)=\frac{\bar{n}(t)}{{{n}_{0}}}. \end{array}\end{eqnarray} \tag{ 179 }$$

In [82], the analysis reveals that for $\beta >0.13$ the probability distribution $p\left(\hat{n}\right)$ at time $t=500\min$ is bimodal with peaks at zero and close to one (figure 9 lower panel). To conclude, simulating the Master equations reveals a continuum of steady state characterized by two peaks, while the mean field approximation predict a bistable distribution, where the dynamics can fall into one of two attractors. More stochastic analysis is expected to reveal in the future how gene expression regulate development, boundary between brain regions [83] or diseases.