Timing and Shape of Stochastic Autocatalytic Burst Formation

TL;DR

This paper investigates how stochastic noise influences the timing and shape of burst formation in autocatalytic systems near a saddle-node bifurcation, revealing pattern nucleation driven by microscopic fluctuations.

Contribution

It applies Large Deviation Theory to characterize burst shape and timing in spatially extended stochastic autocatalytic systems near bifurcation points.

Findings

Burst shape and timing scale with noise amplitude

Spatial nucleation of peaks occurs randomly across the domain

Microscopic noise induces macroscopic pattern formation

Abstract

Chemical, physical and ecological systems passing through a saddle-node bifurcation will, momentarily, find themselves balanced at a semi-stable steady state. If perturbed by noise, such systems will escape from the zero-steady state, with escape time sensitive to noise. When the model is extended to include space, this leads to different points in space "escaping from zero" at different times, and uniform initial conditions nucleate into sharp peaks spread randomly across a nearly uniform background, a phenomena closely resembling nucleation during phase transition. We use Large Deviation Theory to determine burst shape and temporal scaling with respect to noise amplitude. These results give a prototype for a particular form of patternless symmetry breaking in the vicinity of a stability boundary, and demonstrates how microscopic noise can lead to macroscopic effects in such a region.

Tables

Table 1. Table 1: Comparison of action functionals and approximate escape times.

Equation	$m\left(T\right)$	$𝔼\left(T\right)\sim$
$u_t=u^2+ϵ\xi$	$O\left(T^{-3}\right)$	${ϵ}^{-2/3}$
$u_t=u^2+{\gamma }^2+ϵ\xi$	0	${\gamma }^{-1}$
$u_t=u^2-{\gamma }^2+ϵ\xi$	$16{\gamma }^3/3$	$\exp \left[{ϵ}^{-2}16{\gamma }^3/3\right]$
$u_t=\Delta u+u+ϵ\xi$	$U_{f}O\left(e^{-2T}\right)$	$\log \left(ϵ\right)$
$u_t=\Delta u+u^2+ϵ\xi$	$O\left(T^{-5/2}\right)$	${ϵ}^{-4/5}$

Equations

$$u_t=u_{xx}+u^2+ϵ\xi ,$$

$$u_t=D\Delta u+f(u,v)-\delta u$$

$$\text{where }f(u,v)=v\left(k_0+\frac{\gamma u^2}{K^2+u^2}\right),$$

$$\dot{u}=u^2+ϵ\xi .$$

$${\bar{S}}_T(\xi )={\int }_0^T{ϵ}^2{\xi }^{2}dt.$$

$$S_T(\varphi )={\int }_0^T(\dot{\varphi }-{\varphi }^2{)}^{2}dt$$

$$\underset{ϵ\to 0}{\lim }{ϵ}^2\ln \mathbb{E}(\tau )=\underset{T,\varphi }{\min }{\bar{S}}_T(\varphi ),$$

$$\underset{ϵ\to 0}{\lim }{ϵ}^2\ln P(\tau \le T)=-\underset{\varphi }{\min }{\bar{S}}_T(\varphi )$$

$$0=\frac{d{\bar{S}}_T(\varphi )}{d\varphi }=-2\ddot{\varphi }+4{\varphi }^3,$$

$${\int }_{U_i}^{\varphi (t)}\frac{1}{\sqrt{{\psi }^4+C}}d\psi =T-0.$$

$$C\approx {\left(\frac{4\Gamma (5/4{)}^2}{\sqrt{\pi }(T+1/U_f)}\right)}^4.$$

$$S_T(\varphi )$$

$$\approx \frac{C}{3U_f}-0+\frac{1}{3}CT\text{ , for large Uf}.$$

$$\underset{ϵ\to 0}{\lim }{ϵ}^2\ln [P(\tau \le T)]$$

$$\Rightarrow P(\tau \le T)$$

$$\mathbb{E}(\tau )$$

$$=O({ϵ}^{-2/3})$$

$$\dot{u}=u^2+{\gamma }^2+ϵ\xi$$

$$\dot{u}=u^2-{\gamma }^2+ϵ\xi$$

$$u_t=u_{xx}+u^2+ϵ\xi .$$

$$u_t=u_{xx}+u+ϵ\xi .$$

$$S_T(u)={{\int }_0}^T{{\int }_{-\infty }}^{\infty }(u_t-u_{xx}-u{)}^{2}dxdt,$$

$$\begin{array}{r}{u}_t=u_{xx}+u+ϵ\xi ,\\ u(x,0)=0,\\ -{\xi }_t={\xi }_{xx}+\xi ,\\ \xi (x,T)=\alpha \delta (x).\end{array}$$

$$\xi (x,t)$$

$$u(x,t)$$

$$\begin{array}{rl}\text{Minimizing }{S}_T(u)& ={{\int }_0}^T{{\int }_{-\infty }}^{\infty }(u_t-u_{xx}-u^2{)}^{2}dxdt\end{array}$$

$$\begin{array}{r}{u}_t=u_{xx}+u^2+ϵ\xi ,\\ u(x,0)=0,\\ -{\xi }_t={\xi }_{xx}+2u\xi ,\\ \xi (x,T)=\alpha \delta (x).\end{array}$$

$$-{\hat{\xi }}_t={\hat{\xi }}_{xx}+2u_0\hat{\xi }$$

$$\mathbb{E}||u|{|}^2=\sum ϵ(1-e^{-2(1+|\mathbf{k}{|}^2)t})/2(1+|\mathbf{k}{|}^2).$$

$$U_f=\frac{\alpha ϵ}{2}{\int }_0^{2T}{e}^{\gamma }N(0,\gamma )d\gamma ,$$

$$u_i=0.9u_{i-1}+ϵ{\xi }_i,$$

$$u_0=0.$$

$$P(X)={\int }_{u\in X}{p}_u(\mathbf{u})d\mathbf{u}={\int }_{\xi \in X}{p}_{\xi }(\mathbf{\xi })d\mathbf{\xi }.$$

Full text

1 Introduction

The universe we live in does not consist of a perfectly uniform cloud of hydrogen gas, nor is the earth a perfect sphere. Clouds condense into droplets of rain [McGraw and Liu(2003)], a single zygote divides and differentiates into a multitude of different cell types, a flipped coin lands either heads, or tails. All around us, symmetry is broken every day.

Depending on the system of interest, symmetry can be broken in either a patterned or unpatterned manner. Pattern formation is ubiquitous in nature, from a Zebra’s strips [Turing(1952)] to Tiger bush [Lefever and Lejeune(1997)]. Patterned symmetry breaking in reaction diffusion systems is well understood, and was first studied by Alan Turing in the 1950s [Turing(1952)].

In nature we also observe patternless symmetry breaking, in for example the nucleation of crystals[Erdemir et al.(2009)Erdemir, Lee, and Myerson] or rain drops[McGraw and Liu(2003)]. This patternless symmetry breaking is less studied in a mathematical context, and is the focus of this article.

Our present work was motivated by intriguing experimental results in a recent study of the Min cell division system in E.coli [Vecchiarelli et al.(2016)Vecchiarelli, Li, Mizuuchi, Hwang, Seol, Neuman, and Mizuuchi]. In their study, Vecchiarelli et al. observed that Min proteins demonstrate either patterned or unpatterned symmetry breaking, depending on the concentration conditions. For high concentrations of the relevant Min proteins, spirals and linear waves are observed, with each of these patterns possessing their own well studied mathematical frameworks [Murray(2008)]. For lower concentrations, Vecchiarelli et al’s experiment also demonstrates “burst” behavior, whereby the local concentration of Min proteins increases for spatially localized patches of membrane, and then reduces again shortly thereafter (see figure 1). These bursts appear to arise from a relatively homogeneous background, and in contrast to the ordered high concentration patterns described above, they are positioned randomly across the membrane, and appear governed by low level particle scale noise.

In this paper, we demonstrate how dynamical systems in the vicinity of a semi-stable steady state are sensitive to low level white noise, giving rise to similar burst formation behavior. We consider a “canonical” stochastic PDE model of such a system

[TABLE]

and use Large Deviation Theory to predict both the time of burst formation, and the shape of the burst formed, in the limit of low amplitude noise.

We consider a particular reaction-diffusion system from the literature in Sect. 2, and use simulations to demonstrate the occurrence of burst formation in appropriate parameter regimes. In Sec. 3, we use Large Deviation Theory to predict the time taken for the stochastic ODE $u_{t}=u^{2}+\epsilon\xi$ to evolve from $u=0$ to some fixed value $U_{f}$, and determine the most probable path for this occurrence. Sect. 4 extends this work to systems with one spatial dimension via the addition of a diffusion term (as is 1). We use a mixture of numeric and analytic techniques to determine both the shape and time taken for bursts to form. In Sec. 5, we compare our analytic results both to one another, and to the corresponding simulations, and in Sec. 6 we discuss the conceptual difficulties preventing our analysis from being extended to higher dimensions. Finally, in Sec. 7 we summarize the results, and discuss limitations, physical relevance, and future directions for research. The Appendices contain both an introduction to Large Deviation Theory for those unfamiliar with the topic, along with more detailed calculations excluded from the main discussion.

2 Computational Exploration

Before delving into analytical results, we take a brief computational detour to demonstrate the general behaviour of the class of neutrally stable systems that form the focus of this paper.

Consider the following equation due to Mori et al. [Mori et al.(2008)Mori, Jilkine, and Edelstein-Keshet]

[TABLE]

over the unit square in $\mathbb{R}^{2}$ with periodic boundary conditions. For the sake of these simulations, we use $\gamma=5$, $K=2$ $\delta=0.5$ and diffusion coefficient $D=10^{-4}$.

Equation (2) was designed to model a system of cell signaling proteins called Rho GTPases. The model was originally used to demonstrate wave pinning, a behavior reminiscent of the “cell polarization” behavior associated with Rho proteins.

In this equation $D$, $k_{0}$, $\gamma$, $K$ and $\delta$ are reaction and diffusion parameters, $u$ is the concentration of membrane bound Rho protein and $v$ is the concentration of Rho diffusing quickly through the bulk of the cell. As $v$ is varied, the above system gains and loses steady states via a pair of saddle-node bifurcations (see figure 3). Let $u_{0}$ and $v_{0}$ denote the values of $u$ and $v$ at the second of these saddle-node bifurcations, in which the lower stable steady state is annihilated. In this region $v_{0}\approx 0.5$.

Simulating with fixed $v=v_{0}$, $u(0,{\bf x})=u_{0}$ and additive noise of the form $\epsilon\xi$, leads to a long nearly homogeneous period where $u({\bf x},t)\approx u_{0}$ for all ${\bf x}$, followed by the abrupt formation of several “spikes”, which then saturate and spread across the domain. (See figure 3 A)

By contrast, when we select some fixed $v$ greater than $v_{0}$ the system remains approximately homogeneous as it tend towards this $u$-saturated state. For $v$ less than $v_{0}$, the system instead approaches some homogeneous low $u$ state. In either of these latter two cases, no significant patterning or heterogeneity is observed. (See figure 3 B & C)

Similar results may be observed for a variety of spatially distributed systems in the vicinity of saddle-node bifurcations. In order to consider this behavior in more generality, we consider a variety of “cannonical” equations, each demonstrating different aspects of the behavior of interest in the vicinity of a generic saddle-node.

3 Escape time from steady states; spaceless case

Let us begin our discussion by considering first the spaceless stochastic differential equation

[TABLE]

Here $u$ is our stochastic process of interest, $\epsilon$ is some small noise amplitude parameter, and $\xi$ is assumed to be Gaussian white noise. We assume noise as defined by Walsh [Walsh(1986)], such that the integral over any time interval $\Omega$ is given by $\int_{\Omega}\xi dt=N(0,|\Omega|)$, and integrals over non-intersecting time intervals are independent. The extension to spatio-temporal white noise (also discussed by Walsh) makes use of higher dimensional integration. Equation 4 provides a canonical example of a spaceless noisy system balanced at a saddle node, and provides a toy model on which to build our understanding before moving to the spatial case.

Deterministically, trajectories starting at $u=0$ remain stationary, and never reach positive $U_{f}$. However, any noise added to the system is enough to disturb this equilibrium and allow $u$ to escape. We would like to determine the amount of time it takes our stochastic process $u$ to travel from the semi-stable steady state $u=U_{i}=0$ to some positive constant $U_{f}$. The appropriate tool for this task is Large Deviation Theory.

Large Deviation Theory (LDT) is a theory used to study unlikely noise driven events in stochastic systems (often, but not always SDEs). For the interested reader, we give a brief introduction to LDT in Appendix A. For those who wish for more detail, we recommend Freidln and Wentzell’s “Random Perturbations of Dynamical Systems” [Freidlin and Wentzell(2012)], and Rassoul-agha and Seppalainen’s “A Course on Large Deviations With an Introduction to Gibbs Measures”[Rassoul-agha and Seppalainen(2015)]. Here we present only the key results of LDT - those concepts and theorems critical to our present interest.

The central premise of LDT is that when unlikely events occur, they are overwhelmingly likely to occur via a path “close to” the most probable path. The most probable path, $\phi(t)$, is found by minimizing the “normalized action functional”

[TABLE]

For any given SPDE, $\bar{S}_{T}(\xi)$ can be re-written as a function of $\phi$. For example, when studying eq. 4 we can rearrange to isolate $\epsilon\xi$ and find

[TABLE]

Minimizing $S_{T}(\phi)$ allows us to determine the most probably path for a given event, and having identified this part we are able to approximate the probability of an event occurring by a given time using either

[TABLE]

to estimate the expected time $\mathbb{E}(\tau)$ or

[TABLE]

to estiamte the probability distribution.

These formula are given as theorem 4.1 and theorem 1.2 (respectively) in chapter 4 of Freidln and Wentzell[Freidlin and Wentzell(2012)].

We will make use of both of these theorems throughout this paper, however, because $\min_{T,\phi}S_{T}(\phi)=0$ for many of the examples discussed, eq. 8 will prove to be the more helpful of the two. This is because we are in the slightly unusual position of studying escape from an energy plateau, as opposed to an energy well, as might be more usually studied.

Suppose we wish to study the probability of $u$ passing from $U_{i}=0$ to $u\geq U_{f}$ by time $t=T$, that is to say the “escape from zero” problem for equation 4. In order to use either of the above formula, we must find the minimum of ${S}_{T}(\phi)$. By the calculus of variations, any minimizer of ${S}_{T}(\phi)$ must satisfy

[TABLE]

and hence

[TABLE]

For $C=0$ we recover the deterministic solution $\phi(t)=1/(U_{i}^{-1}-t)$. This solution indicates infinite travel time to $U_{f}$ in the case where $U_{i}=0$ (our case of interest). For $C\neq 0$, the integral can be solved using wolfram|alpha [Wolfram|Alpha(2018a)], determining $t$ in terms of $\phi$. Given that the function itself provides limited illumination, and is defined in terms of the ellipticF function (itself defined in terms of an integral), we will refer to solutions of 10 simple as $\tau(\phi)$, the inverse of $\phi(t)$.

In order to match our assumed boundary conditions, we must pick $C$ such that $\tau(U_{f})=T$. This is not generically an easy problem to solve exactly, but for large $U_{f}$ it can be well approximated. In order to determine $C$, in this limit, we first note $\tau_{\infty}=\lim_{\phi\rightarrow\infty}\tau(\phi)=\frac{4\Gamma(5/4)^{2}}{\sqrt{\pi}C^{1/4}}$ (as given by wolfram|alpha [Wolfram|Alpha(2018b)]). Here $\Gamma$ is the Gamma function [Artin(2015)], and takes the value $\Gamma(5/4)\approx 0.9064$. The deterministic time taken to get from $\phi=U_{f}$ to $\phi=\infty$ is $1/U_{f}$, thus $\tau_{\infty}\approx T+1/U_{f}$. Rearranging gives

[TABLE]

Numerical experimentation indicates that, for this choice of $C$, $|\tau(U_{f})-T|/T<10^{-6}$ whenever $T,U_{f}\geq 4$, indicating that our approximation of $C$ is very good.

Having determined $C$, we are able to plot $t=\tau(\phi)$ for arbitrary $U_{f}$ and $T$ (see figure 4, Left).

In order to invoke either eq. 7 or eq. 8, we must must determine $S_{T}(\phi)$ along this action minimizing path. Starting from the definition of $S_{T}(\phi)$ (equation 6), it is possible to show that in our case

[TABLE]

Details are given in appendix B. Remembering from equation (11) that $C=O(T^{-4})$ we find ${S}_{T}(\phi)=O(T^{-3})$.

This analytic result agrees with numeric results found by approximating ${S}_{T}(\phi)$ directly for a variety of different $T$ values.

Because ${S}_{T}(\phi)\rightarrow 0$ as $T\rightarrow\infty$, we see that eq. 7 gives no information. This is a reflection of the fact that $\phi$ need not overcome any energy barrier in its path from [math] to $U_{f}$, but merely an energy plateau.

We can however apply equation 8 to find

[TABLE]

for small $\epsilon$.

[TABLE]

Hence we find that time taken for a solution to $u_{t}=u^{2}+\epsilon\xi$ to escape from zero and approach infinity is predicted to scale like $O(\epsilon^{-2/3})$. The majority of this time is spent close to $u=0$, and the escape time behavior of the system is dominated by the behavior of the system in this region; the exact value of $U_{f}$ has negligible impact, so long as $U_{f}$ is selected such that $U_{f}\gg O(\epsilon^{2/3})$. Simulation of equation 4 for a variety of $\epsilon$ agrees with these asymptotic results (see figure 4, Right).

3.1 Systems near a saddle-node bifurcation

The above work describes the behaviour of the system precisely balanced at a saddle-node bifurcation. For the sake of completeness, we might also consider systems in the vicinity of such a bifurcation.

The system

[TABLE]

is a canonical example of a system which, for small $\gamma$, has just lost its steady states. In the deterministic limit $\epsilon\rightarrow 0$, this system admits solutions of the form $T=\left[\tan^{-1}(u/\gamma)/\gamma\right]_{U_{i}}^{U_{f}}$, which gives travel times of order $O(\gamma^{-1})$.

The system

[TABLE]

is a canonical example of a system with a stable/unstable pair of steady states. Assuming in this case that $U_{i}=-\gamma$ (the stable steady state), and $U_{f}>\gamma$, it is possible to show that the action minimizing path obeys $\dot{\phi}^{2}=(\phi^{2}-\gamma^{2})^{2}+C^{2}$. In the limit of large $T$, we can select $C\rightarrow 0$, and find $\bar{S}_{T}(\phi)\rightarrow 16\gamma^{3}/3$.

By eq. 7 we thus have $\mathbb{E}(\tau)\sim\exp[\epsilon^{-2}16\gamma^{3}/3]$.

In addition to the above static results, the behavior of systems passing through saddle node bifurcations (either with or without noise), has also been studied. A review of this “delayed bifurcation theory” is given by Christian Kuehn [Kuehn(2011)]. For a more detailed introduction, see Berglund and Gentz [Berglund and Gentz(2006)]. We now move on to consider the previously unexplored spatially distributed problem.

4 The Spatially Distributed Problem

Now that we have built up our understanding using the spaceless model, let us turn our attention to the spatially distributed model. Our primary interest in the spatially distributed case is considering the behavior of systems at a saddle-node bifurcation. Spatial systems of this form can be represented via the canonical equation

[TABLE]

Here, both $u$ and $\xi$ are functions of $x$ and $t$. We assume that $u$ starts at the steady state; $u(x,0)=0$. Here $\xi$ represents our white noise term, although it may be thought of as a forcing function, with an associated energy $\int\xi^{2}dxdt$.

We wish to describe the behavior of the system as it ‘escapes’ from the steady state at zero. We assume (without loss of generality) that $u$ takes its maximum at $x=0$, and ask “for a given time $T$, what is the most probable path such that $u(0,T)=U_{f}$, where $u(x,t)<U_{f}$ for all $t<T$?”

For the sake of comparison, we will also consider escape from a linear unstable stead state,

[TABLE]

While less mathematically interesting, this linear escape problem gives us an analytically tractable case to explore, along with something to compare our non-linear results to. We begin with this simpler case.

4.1 Linearly Unstable Case

In what follows we present only the most important results, and the conceptually important steps leading to these results. Details can be found in Appendix C. In order to understand spike formation in equation (21), we must minimize the associated “action functional”:

[TABLE]

subject to the conditions $u(x,0)=0$ and $u(0,T)=U_{f}$. As per standard functional analysis techniques, this minima can be found when the functional derivative $\frac{dS_{T}(u)}{du}(x,t)=0$. This condition can be written as a single differential equation with double the number of derivatives with respect to each variable, but is more conveniently written as a coupled pair of PDEs:

[TABLE]

This system admits the explicit solution:

[TABLE]

where $N(x,\gamma)=\frac{1}{\sqrt{4\pi\gamma}}\exp\left[\frac{-x^{2}}{4\gamma}\right],$ and $\alpha=\frac{4U_{f}}{\epsilon erfi(\sqrt{2T})}$.

The associated action functional can be shown to be $S_{T}(u)=4U_{f}^{2}/erfi(2T)=U_{f}^{2}O(e^{-2T})$.

Further analysis and discussion of these results we postpone until section 5.1.

4.2 Semi-stable Case

We now consider the non-linear case, (equation 20). As previously, in order to determine the dynamics of spike formation, we concern ourselves primarily with determining the most probable path to spike formation. We present here only the conceptually important milestones in our calculations. Details are similar to those given in Appendix C.

[TABLE]

subject to the constraints $u(x,0)=0$ and $u(0,T)=U_{f}$, we find that the minimizer $u(x,t)$ satisfies

[TABLE]

Unlike the linear case, where $\xi$ could be solved independently of $u$, and then used to determine $u$ explicitly, here no such analytic solution is available. $\xi$ and $u$ are inextricably coupled. Solving 27 numerically in either direction in time requires that we solve the ill posed backwards heat equation, either for $u$ or $\xi$, depending on which direction we solve in. Creating a full spatial-temporal mesh and solving for both $\xi$ and $u$ simultaneously is possible, but quickly becomes computationally expensive for finer meshes. We can, however, solve iteratively. This is done by initially assuming $u_{0}(x,t)=0$, and solving

[TABLE]

with the boundary condition $\hat{\xi}(x,T)=\delta(x)$, to determine $\hat{\xi}_{1}$. We use $\alpha\hat{\xi}_{1}$ as a forcing term to find $u_{1}$, and $u_{1}$ to determine $\hat{\xi}_{2}$ and so on, until $|S(u_{i})-S(u_{i+1})|<0.0005|S(u_{i})+S(u_{i+1})|$. (See figure 5). Using this iteration scheme, we are able to determine $u,\xi$ and $S_{T}(u)$ for the action minimizing $u$ profile. The value of $\alpha$ in each iteration step is found using a bisection method in which $\alpha$ is reduced if $u(0,T)>U_{f}$ and increased if $u(0,T)<U_{f}$.

5 Comparisons of results

5.1 Comparison Between Results for Linear and Non-Linear Case

Now that we have determined the action minimizing profile in both the linear and non-linear case (equations (21) and (20), respectively), let us now compare these two models.

To begin, we would like to determine the time till burst formation in each model. For the sake of notation, let us define $m(T)=\min_{\phi}S_{T}(\phi)$. While under normal circumstances we would find the expected escape time by invoking equation 7, in all cases studied here, $m(T)\rightarrow 0$ as $T\rightarrow\infty$. Hence equation 7 provides us with no information, and we instead rely on equation 8.

Regardless of the particular system under study, when $m(T)\gg\epsilon^{2}$ equation 8 implies $-\ln P(\tau\leq T)\gg 1$ and hence $P(\tau\leq T)\approx 0$. This in turn implies that the probability density of burst formation for such values of $T$ is low. Similarly, when $m(T)\ll\epsilon^{2}$ we have $-\ln P(\tau\leq T)\ll 1$ and hence $P(\tau\leq T)\approx 1$, and once again the probability of escape is negligible, as escape has almost certainly already occurred. Escape can only have non-negligible probability density when $\ln P(\tau\leq T)$ is neither too big nor too small, namely when $m(T)$ and $\epsilon^{2}$ are similar orders of magnitude.

This implies that for the linear model escape is predicted when $\epsilon^{2}\approx\bar{S}_{T}=U_{f}^{2}O(e^{-2T})$, or equivalently when $T\sim-\log(\epsilon)$. The time of escape changes little, even as $\epsilon$ is varied over several orders of magnitude. In the nonlinear case we observe that $m(T)=O(T^{-5/2})$ as $T\rightarrow\infty$ (see figure 6). By the above arguments, we predict $T\sim\epsilon^{-4/5}$; the time till burst formation is sensitive to the amplitude of the noise driving the system. Table 1 summarises these timing results, and compares to the previously discussed spaceless models.

Far more interesting than the difference in escape times for the linear and quadratic case is the substantial difference in the profile of our action-minimizing $u$. In the linear case, $u(x,t)$ was found explicitly (equation 25). In this case $u(x,T)$ is well approximated by a corresponding normal distribution $N(x,2T)$, and becomes wider over time (see figure 7, upper). In contrast, in the quadratic case, the final width of $u$ increases at most very slowly, and appears to approach a limiting distribution as $T\rightarrow\infty$ (see figure 7, lower).

As a final point of comparison, we consider the reaction of the two systems to changes in $U_{f}$. As might be expected, the shape of the action-minimising $u$ profile is unaffected by changes in $U_{f}$ for the linear system. Doubling $U_{f}$ for a given $T$ doubles both $\xi$ and $u$, but does not change the shape of either. In contrast, in the non-linear case, the profile of $u(x,t)$ is sensitive to changes in $U_{f}$, with larger $U_{f}$ leading to a sharper “spike” solution, and smaller $U_{f}$ leading to soft “bumps” (see figure 8).

5.2 Comparison of Analytic Results to Direct Simulations

As a final check on the above analytic results, we compare to simulations of the corresponding systems.

In the linear case (equation 21), we find that the time until $\max|u(x,t)|=U_{f}=1$ scales like $-\log(\epsilon)$, as predicted by theory (see figure 9, left). In the non-linear case (equation 20), burst time is reasonably approximated by $\epsilon^{-4/5}$ as predicted by theory, but the agreement is weaker than might be hoped (see figure 9, right). It is unclear if this disagreement stems from the limitations of our simulations, or the numeric-analytic arguments proposed earlier.

Comparison of shape in the linear case finds reasonable agreement in the vicinity of the primary burst itself, however because a large number of bursts form simultaneously, the final profile of $u$ ends up being an overlapping combination of many different such bursts (see figure 10,left).

In contrast, comparison between the predicted and observed shape of bursts shows strong agreement in the non-linear case, even for relatively large values of $\epsilon$ (see figure 10,right). Bursts spend a greater portion of their development time close to zero in the non-linear case, and hence, although many bursts may begin growing at roughly the same time, the second and third place bursts will have only a minor effect on the overall profile of $u(x,T)$.

The code for all simulations can be found in the supplementary materials.

6 Difficulties in 2D, and higher dimensions

Given that the universe we occupy is not one dimensional, it would be beneficial to extend the above results to higher dimensions, preferably three dimensions, although for certain contexts, such as the binding of proteins to a cell membrane[Vecchiarelli et al.(2016)Vecchiarelli, Li, Mizuuchi, Hwang, Seol, Neuman, and Mizuuchi], two dimensions would suffice.

Unfortunately, the properties of white noise preclude this avenue of investigation; as noted by Ryser et al [Ryser et al.(2012)Ryser, Nigam, and Tupper], equations of the form $u_{t}=\Delta u+\xi$ do not remain well posed in dimension two or higher. The smoothing effects of diffusion are insufficient to restrain the irregularity produced by space-time white noise.

This result can be seen through a variety of lenses, depending on ones outlook. Ryser et al, describe $u$ in terms of a sum of Fourier modes, and show

[TABLE]

In one dimension, when ${\bf k}\in\mathbb{Z}$ this sum converges, however for two or more dimensions, ${\bf k}\in\mathbb{Z}^{d}$ and the sum is unbounded for all $t>0$. This indicates that $u$ is not well defined. Simulations in 2d, such as those presented in section 5.1 function for fixed $dx$, but invariably break down as $dx\rightarrow 0$.

In terms of the methodology used in this paper, the degeneracy of higher dimensions manifests as a break down of various integrals. Issues arise regardless of what equation is used, but the problem is most easily illustrated by considering the linear case. Consider the solution $u(x,t)$ of eq. 21 given in eq. 25. The parameter $\alpha$ of eq. 25 is determined by considering the boundary condition $u(0,T)=U_{f}$,

[TABLE]

where $\gamma=T+t-2\tau$.

In 1D $N({\bf x},\gamma)=(4\pi\gamma)^{-1/2}\exp\left[-x^{2}/4\gamma\right]$, and hence $U_{f}=\frac{\alpha\epsilon}{4\sqrt{\pi}}\int e^{\gamma}\gamma^{-1/2}e^{-0/\gamma}d\gamma$. Because the integral converges to a finite number, we can rearrange to determine $\alpha$ and go on to find $S(u)$ explicitly for any $U_{f}$.

By contrast, in 2D we have $N({\bf x},\gamma)=(4\pi\gamma)^{-1}\exp\left[-\|{\bf x}\|^{2}/4\gamma\right]$, and hence, $U_{f}=O(\alpha\epsilon)\int e^{\gamma}\gamma^{-1}e^{-0/\gamma}dx$. The integral is unbounded, and thus we are able to select $\alpha$ arbitrarily small. This in turn allows both $T$ and $S(u)$ to be selected arbitrarily close to zero, regardless of $U_{f}$. Physically, we are picking $\xi$ to be a spike that is sufficiently narrow so as to ensure $S(u)$ arbitrarily small, yet also sufficiently high so as to force $u$ from zero to $U_{f}$ in time $T$. This sleight of hand is impossible in one dimensional systems, but for all higher dimensions it is the “action minimizing” strategy.

While here we have described the difficulties of integration in the linear case, the problems described extends to all equations that include noise terms in 2D, and become progressively worse for higher dimensions. They are a symptom, rather than the root cause of the problems of noise in higher dimensions.

7 Conclusions, Discussion and Future work

In this article we have presented a number of prototypical models exploring the behavior of systems starting at unstable and semi-stable steady states.

We have demonstrated both analytically and through simulations that under the influence of spatio-temporal white noise, systems of the form $u_{t}=u_{xx}+u^{2}+\epsilon\xi$ develop well defined “spikes” in finite time. This phenomena can be seen as delayed bifurcation [Kuehn(2011)] in a spatially distributed system, and also as a compliment to non-linear blow up phenomena, in which we consider instead the system’s “escape from zero”, as opposed to the more typically studied “approach to infinity” [Vázquez and Galaktionov(2002)].

The symmetry breaking observed here does not rely on multiple chemical species, nor on contrasts in diffusion rate, as might be expected for the more classical Turing type pattern formation.

We were able to determine both the “energy cost” of spike formation, along with the shape of the most probable spike, and showed that in the limit of large time, spike formation requires infinitesimally small energy. The sensitivity of the system in the vicinity of its semi-stable steady state allows noise of amplitude $\epsilon\ll 1$ to have macroscopic effects in O($\epsilon^{-4/5}$) time. We also explored similar results both in the spaceless case, and in the case where our $u=0$ state is linearly unstable as opposed to semi-stable.

In terms of physical relevance, the systems studied can be best thought of as a prototype for chemical reaction-diffusion systems in the vicinity of a saddle-node bifurcation. Our study demonstrates the breakdown of a homogeneous state, and can be seen as a first step towards understanding systems such as the Min system studied by Vecchiarelli et al [Vecchiarelli et al.(2016)Vecchiarelli, Li, Mizuuchi, Hwang, Seol, Neuman, and Mizuuchi], in which changes in the bulk concentration of a protein pushes the system through a bifurcation boundary and leads to the formation of membrane bound “bursts”.

There are a number of further questions which must be answered before the work here can be applied to any experimental context. First and foremost, the work here considers a system perfectly balanced at the saddle-node; $u_{t}=u^{2}$. No physical system however is ever so perfectly balanced, and so determining the robustness of these results for a system that is passing through such a saddle node bifurcation is critical to our understanding of the relevance of these results. Further in order for the results here to prove useful for experimental scientists, time must be invested in investigating the mapping between experimental observations, and the associated reaction and diffusion parameters implied; for example, does the relationship between spike width and height give a reliable signature that can be used to determine system parameters?

Finally, our work here alludes (in passing) to potentially deeper questions in chemistry; namely the difficulties in standard representations of noise when in higher dimensions. Such questions must be answered, or at the very least sidestepped, before the phenomena observed here can be sensibly applied to the multi dimensional systems ubiquitous in the real world.

8 Code

The Code used in this project is available on github at alastair-JL/StochasticBurst.

9 Acknowledgements

We wish to acknowledge Anthony Vecchiarelli, for stimulating discussion on the topic of Min proteins, and sharing his data. This research was funded by the four year fellowship from the university of British Columbia, and NSERC.

Appendix A Introduction to Large Deviation Theory

Here we give a very brief introduction to the principle ideas and techniques used in Large Deviation Theory (LDT), and used throughout this paper. In this appendix we present LDT as it applies to discrete time Stochastic processes. The continuous time case is conceptually similarly. Our goal here is to build intuition rather than mathematical rigor.

As a concrete example with which to frame our discussion, consider the Discrete Orenstien-Ulembeck process[Larralde(2004)]

[TABLE]

Here $\xi_{i}$ is assumed to be a normal random variable with mean zero and variance one, such that each $\xi_{i}$ is independent. We assume $0<\epsilon\ll 1$.

Suppose we wish to know the first time that $u_{i}>1$. At first glance, given $O(\epsilon)$ noise, and the decay term of our OU process, $u_{i}>1$ seems unlikely to occur. That said, as we take $i\rightarrow\infty$, unlikely, even exceedingly unlikely events should occur eventually. LDT concerns itself with such questions as “how long will it take for X to occur?” and “given that X occurs, what path is the system most likely to take in order to get there?”. In our particular case X is “$u_{i}>1$”.

The total probability of our event $X$ is equal to the integral over the probability density of all paths leading to that event. This integral can be formulated either in terms of the paths of the stochastic process ${\bf u}$, or in terms of the underlying noise $\bm{\xi}$.

[TABLE]

The probability density function for a noise vector $\bm{\xi}$ of length $N$ is

[TABLE]

Unfortunately, even with this well defined probability density, the boundary of the integral $\int_{\xi\in X}p_{\xi}({\bm{\xi}})d\bf{\bm{\xi}}$ is generically complicated enough so as to prevent us from evaluating $P(X)$ directly (the OU process being a notable exception).

In order to avoid this complex integration boundary, it is useful to transform our integral back into a form dependent on ${\bf u}$. To achieve this we rearrange our recurrence relation (equation 29), and find

[TABLE]

Using this one-to-one correspondence between ${\bf u}$ and ${\bm{\xi}}$ along with the change of random variables formula, we find

[TABLE]

Here $S({\bf u})$ is said to be the “normalized action functional” of our problem. In the particular case discussed here $S({\bf u})=\sum(u_{i}-0.9u_{i-1})^{2}/2$. In general $S({\bf u})$ depends both on the equation governing a stochastic process, and the particular form of the noise generating it. $S({\bf u})$ can be thought of as a measure of the “total improbability” associated with a given path, and is associated with the amount of “energy” that noise must pour into the system in order to cause a particular path to occur. When dealing with continuous systems, $S$ is defined as an integral over the square of noise, rather than a sum.

At this stage, in order to determine the probability of our event, we need to evaluate $\int_{u\in X}e^{-S({\bf u})/\epsilon^{2}}d{\bf u}$. Typically, this integral can not be evaluated exactly, but it can be well approximated via Laplace’s Principle [Laplace(1986), Olivieri and Vares(2005)].

Laplace’s Principle states that:

[TABLE]

where here we minimize $S({\bf u})$ over all ${\bf u}\in X$. Laplace’s principle is based on the idea that for integrals of the from $\int_{u\in X}e^{-S({\bf u})/\epsilon^{2}}d{\bf u}$, the overwhelming majority of all probability mass is concentrated in the vicinity of $\min S({\bf u})$ whenever $\epsilon\ll 1$, (see figure 11). From the point of view of stochastic processes, what Laplace’s Principle is effectively stating is that if an improbable event does occur, the observed path is overwhelmingly likely to be ‘close’ to the most probable path.

Let us return to the original question proposed at the start of this section: how long does it take before the Discrete Orenstien-Ulembeck process defined in equation 29 exceeds one for the first time?

Suppose we wish to determine the probability that $u_{N}\geq 1$, for some particular $N$, assumed to be large. By Laplace’s Principle, we must thus minimize $S({\bf u})=\sum(u_{i}-0.9u_{i-1})^{2}/2$. Because any $u_{N}>1$ will only increase $S({\bf u})$ we can assume $u_{N}=1$. In order to minimize, we require that $u_{N}=1,u_{0}=0$ and $\frac{dS}{du_{i}}=0$ for all $i$. As $N\rightarrow\infty$, solutions to the above can be well approximated by $u_{i}=0.9^{N-i}$. Substituting back into the definition gives $S({\bf u})=0.095$, and hence $P(u_{N}\geq 1)\approx\exp[-0.095/\epsilon^{2}]$. If this is the probability of success for any particular large $N$, then we infer that the expected time until $u_{i}>1$ will scale such that $\mathbb{E}(\tau)=O(e^{0.095/\epsilon^{2}})$.

Freidln and Wentzell[Freidlin and Wentzell(2012)] provide the general formulation for the above two results. They state (chapter 4, theorem 1.2) that for small epsilon the probability of a particular rare event occurring by a given time $P(\tau<T)$ is governed by:

[TABLE]

Here, as is customary, $\phi$ denotes a particular trajectory of $u$ and $S_{T}(\phi)$ indicates that we are minimising over all trajectories $\phi$ such that our rare event occurs by time $T$. This theorem is referred to as equation (8) in the main text.

Appendix B Derivation of analytic expression for $\bar{S}_{T}$ in the spaceless case.

Here we present the derivation of equations 13. We begin by taking the definition of $\bar{S}_{T}$ and expanding. From eq. 9 we can show $\dot{\phi}=\sqrt{\phi^{4}+C}$ and hence:

[TABLE]

The problematic term here is the first integral. Applying integration by parts gives

[TABLE]

By eq. 10 this last term is $2CT$. Substituting this result back in gives.

[TABLE]

hence recovering equation 13.

At no point in this derivation have we used the approximation of $C$ given in eq 11, hence this is an analytic result.

Appendix C The linearly unstable case; full calculation

Here we provide the detailed algebra suppressed in section 4.1. Consider the minimization problem:

[TABLE]

This is equivelent to equation (26) of the main text.

If we consider $\xi$ as a form of forcing function, then the above asks for the lowest energy forcing required to lift $u(0,T)$ to height $U_{f}$. Thinking of $\xi$ as a form of noise we are asking for the most probable noise. For the rest of this appendix we will think of $\xi$ in terms of energy input, although we will keep the white noise formulation in mind. Please note that while both $u$ and $\xi$ are stochastic processes, in a slight abuse of notation we also use the symbols to indicate the optimal path subject to the above constraints.

In order to find our minima, we first re-write $S$ in terms of $u$, and then take the functional derivative:

[TABLE]

[TABLE]

Now using integration by parts on each $\psi$ term, as appropriate, we find

[TABLE]

Because we seek to minimise $\int\int\xi^{2}dxdt$, we can safely demand that $\xi,\xi_{x}\rightarrow 0$ as $x\rightarrow\pm\infty$ and thus $\left[\psi\xi_{x}-\psi_{x}\xi\right]_{-\infty}^{\infty}=0$. We are left with the potentially problematic $\int\left[\psi\xi\right]_{0}^{T}\xi dx$ term, however, since $u$ is fixed at zero for $t=0$ and at $U_{f}$ for $x=0,t=T$ we know $\psi$ (our perturbation to $u$) is zero in these location. Further, we can assume that $\xi=0$ for $t=T,x\neq 0$, as any non-zero forcing at these locations would increase $S$, but have no effect on $u(0,T)$. Hence $\int\left[\psi\xi\right]_{0}^{T}\xi dx=0$, and so

[TABLE]

The above can be stated as a fourth order differential equation purely in terms of $u$, but this provides little illumination. Instead we leave it as a backwards heat equation, ready to be coupled to C.

As argued previously, for our optimal solution $\xi(x,T)=0$ whenever $x\neq 0$, as non-zero values of $\xi$ in these locations increase $S(\xi)$ but have no effect on $u(0,T)$. In order for $\xi$ to have any effect on $u$ it must have non-zero total mass, and since $\int-\xi_{t}dx=\int\xi dx$ we know that it must have non-zero total mass at time $T$. Hence $\xi(x,T)=\alpha\delta(x)$ for some $\alpha$.

Combining the above, the minimal noise carrying $u$ from [math] to $u(0,T)=U_{f}$ must solve:

[TABLE]

This recovers equation (27) of the main text.

Solving directly gives

[TABLE]

Where

[TABLE]

The noise scaling $\alpha$ can be found using our boundary condition:

[TABLE]

In order to find the normalized action functional we find:

[TABLE]

Escape is predicted to occur when $\epsilon^{-1}=O(e^{2T})$. For all but the smallest $\epsilon$ values, this can be considered to take place in $O(1)$ time.

Conceptually similar calculations can be used to get from equation (26) to (27) in the non-linear case.