Brownian bridges for stochastic chemical processes

Journal club presentation by Kyunghoon Han

Theoretical Chemical Physics Group, University of Luxembourg

Some facts on the paper

Title : Brownian bridges for stochastic chemical processes--An approximation method based on the asymptotic behavior of the backward Fokker-Planck equation

Authors : Shiyan Wang, Anirudh Venkatesh, Doraiswami Ramkrishna, Vivek Narsimhan of Purdue University

Their goal?

A fast approximation method to generate Brownian bridge process without solving the Backward Fokker-Planck (BFP) equation

The game plan

  1. Mathematical background
  2. Problem described by the authors
  3. Results and discussion
  4. Relevance to my work?

Brownian motion

A Brownian motion is described by what is called a Wiener process, $W_t$, parametrized by a variable $t$, which is defined as:
  1. $W_0 = 0$
  2. $W_t-W_s \sim \mathcal{N}(0,t-s)$
  3. Each increment of $W_t$ is independent and is almost surely continuous

Brownian Bridge?

Suppose $(\mathcal{M},g)$ is a Riemannian manifold with $\dim\mathcal{M} = n$. If $X_t(s)$ is a Brownian motion defined on $(\mathcal{M},g)$, parametrized by a variable $t$ such that $X(0)=x\in\mathcal{M}$. The probability of finding $X_s(x)$ in some region $A\subset \mathcal{M}$, then, one may define: $$ \mathbb{P}\left\{X_s(x)\in A\right\} = \int_{\mathcal{M}} 1_{A}(z)p(s,x,z) \text{ vol}(dz) $$
Note that the quantity defined above is essentially computing the expectation of finding the quantity in $A$. The map $1_A$ is an indicator function.
Then the Brownian bridge from $x\in\mathcal{M}$ to $y\in\mathcal{M}$ with lifetime $T\geq 0$, denoted as $X^T(x,y)$, is defined as a Brownian motion $X$ with $X_0=x$ and $X_T=y$ such that \[\begin{aligned} &\mathbb{P}\left\{X_s^T(x,y)\in A\right\} \\ &= \int_{\mathcal{M}} 1_A(z) \frac{p(s,x,z)p(T-s,z,y)}{p(T,x,y)}\text{ vol}(dz) \end{aligned}\]
Essentially, the Brownian bridge considers a stochastic process where the random variable $X_t$ is tied down at $t=0$ and $t=T$.

Simple simulation of a Brownian bridge

5 batches and 50 steps

10 batches and 100 steps

20 batches and 1500 steps

50 batches and 20000 steps

Stochastic D.E. in Itō form

Givens

  1. $s$ - a parameter of the path taken ($0 < s < L$)
  2. $\vec{x}(s)$ - position along the path
  3. drift term $\vec{A}$: the term that forces the position to change in some expected direction
  4. diffusion term $\vec{B}$: the random changes exerted on $\vec{x}$

$$ d\vec{x} = \vec{A}(s,\vec{x}(s))ds +\vec{B}(s,\vec{x}(s)) \cdot d\mathcal{B} $$
Diffusion tensor: $D^{\mu\nu}=\frac{1}{2}B^\mu B^\nu$

Example: Langevin equation

Suppose $\vec{A}=-ㄷ\nabla U$, $ㄷ$ a drift coefficient, for some continuous dimensionless external potential $U(\vec{x})$. Then the Itō's relation becomes: $$ d\vec{x} = ㄷ\cdot \nabla U ds +\vec{B}\sqrt{ds}\cdot \vec{Z} $$ where $\vec{Z} \sim \mathcal{N}^3(0,\vec{1})$.

Brownian bridge SDE?

An SDE that samples the conditional probability $P\left((\vec{s})|\vec{x}(0)=\vec{x}_0,\vec{x}(L)\in \Omega_L\right)$, i.e. a Brownian bridge $\vec{x}^{Br}$, is:
$$ d\vec{x}^{Br} = d\vec{x} + \underbrace{\vec{B}\vec{B}^T \frac{\partial}{\partial \vec{x}}\left(\ln q\right)}_{\text{additional drift}} ds $$ with $\vec{x}(0)=\vec{x}_0$. Where...

$q(\vec{x},s)$ is a hitting probability.
Meaning that it is likely for the original random walk described by $d\vec{x}=\vec{A}ds + \vec{B} \cdot d\mathcal{B}$ to reach the target region, $\Omega_L$, at $s=L$.
Physically speaking, the extra drift term could be described as an entropic drift that drives the path towards the end region.
The physical entropic free-energy barrier for reaching the end-state in $\Omega_L$ is proportional to $$ T\Delta S \sim -\beta \ln (q) $$

Computation of the hitting probability?

Assuming that $q$ is a probability associated to a Markov chain, the backward Fokker-Planck equation is: \[ \mathcal{L}^bq = \frac{\partial q}{\partial s} + \vec{A} \frac{\partial q}{\partial \vec{x}} + \frac{1}{2}\left(\vec{B}\vec{B}^T\right)\frac{\partial^2 q}{\partial \vec{x} \partial \vec{x}} \]
where $q(\vec{x},s)$ is 1 if $\vec{x} \in \Omega_L$ and is zero if $\vec{x}$ could not reach the target region $\Omega_L$.
The backward integration of the above expression from $s=L$ to $s=0$ gives the hitting probability. But a numerical solution of this BFP is infeasible for high dimensional systems.

Path integral approach?

the probability of attaining a specific path of the SDE $d\vec{x}=\vec{A}ds + \vec{B}\cdot d\mathcal{B}$ is expressed discretely as...
$$ p\{\vec{x}\} \sim e^{-\frac{1}{2\Delta s}\sum_{k=0}^{N-1}\left(\Delta \vec{x}_k -\vec{A}\left(\vec{x}_k\right)\right)\cdot \left(\vec{B}\vec{B}^T\right)^{-1} \cdot \left(\Delta \vec{x}_k -\vec{A}\left(\vec{x}_k\right)\right)} $$ where $ \Delta \vec{x}_k = \vec{x}_{k+1}-\vec{x}_k$.
For a Brownian bridge, we alter $\vec{A}\mapsto \vec{A}^\prime = \vec{A}\left(\vec{x}_k\right) + \vec{u}^d $ with $\vec{u}^d =\left(\vec{B}\vec{B}^T\right)\cdot \nabla_{\vec{x}}\ln(q)$. $$ p\{\vec{x}\} \sim e^{-\frac{1}{2\Delta s}\sum_{k=0}^{N-1}\left(\Delta \vec{x}_k -\vec{A}^\prime\left(\vec{x}_k\right)\right)\cdot \left(\vec{B}\vec{B}^T\right)^{-1} \cdot \left(\Delta \vec{x}_k -\vec{A}^\prime\left(\vec{x}_k\right)\right)} $$ the quantity $\vec{u}^d$ is called the exact drift.
In short, one may write the probability function associated to the Brownian bridge as $$ p^{Br}\left\{\vec{x}\right\} \sim 1(\vec{x}(L)\in \Omega_L)p\{\vec{x}\} $$ where $1$ is the indicator function.

Reweighting

Instead of using $q$ as the hitting probability, consider an approximate hitting probability $\psi$.
The object $\psi$ is defined to approximate the drift $\vec{u}^d$ as $\vec{u}^{ap}$ as: $$ d\vec{x} = \left(\vec{A}+\vec{ap}\right)dt +\vec{B}\cdot d\mathcal{B} $$ where...
the approximated drift is $$ \vec{u}^{ap} = \left(\vec{B}\vec{B}^T\right) \cdot \nabla_{\vec{x}}\ln(\psi) $$ which rewrites the expression of the probability of attaining a path under this approximation is...
Itō's formula then gives \[\begin{aligned} \Delta \ln(\psi) &= \frac{\partial \ln(\psi)}{\partial s}\Delta s \\ &+ \frac{\partial \ln(\psi)}{\partial \vec{x}}\Delta\vec{x} \\ &+ \frac{1}{2}(\vec{B}\vec{B}^T):\frac{\partial^2}{\partial \vec{x}\partial\vec{x}}\Delta s \end{aligned}\]
Noting that $e^{\sum \Delta \ln \psi} \sim 1(\vec{x}(L)\in \Omega_L)$ and taking $\Delta s \rightarrow 0$, one gets $$ p^{ap}\{\vec{x}\} \sim p^{Br}\{\vec{x}\}e^{-\int_0^L\mathcal{G}[\ln \psi] ds} $$ where ...
$$ \mathcal{G}[\ln \psi] = \frac{1}{\psi} \mathcal{L}\psi $$

Ornstein-Uhlenbeck

$$ d\vec{x} = -k\vec{x}ds + \sigma d\mathcal{B} $$ for $k$ being the drift constant and $\sigma$ being the volatility of the random fluctuations.

Geometric Brownian Motion

$$ d\vec{x} = k \vec{x}ds + \sigma \vec{x} d\mathcal{B} $$

Result with Ornstein-Uhlenbeck

Normal bridge vs. introduced approximation

Naive sampling vs. introduced sampling