inordinatum

Physics and Mathematics of Disordered Systems

Posts Tagged ‘brownian motion

Stochastic Logistic Growth

with one comment

A classic model for the growth of a population with a finite carrying capacity is logistic growth. It is usually formulated as a differential equation,

\displaystyle \partial_t n(t)=r\, n(t)\left[1-\frac{n(t)}{K}\right]. (1)

Here n(t) is the size of the population at time t, r is the growth rate and K is the carrying capacity.

The dynamics of eq. (1) is deterministic: The initial population n(0)=n_0 > 0 grows (or decays) towards the constant carrying capacity n(t\to\infty) = K, which is a fixed point of eq. (1). This is seen in the solid trajectories in the figure below:

stoch_log_growth

Deterministic vs. stochastic logistic growth. Solid lines: deterministic model in eq. (1). Transparent lines: Realizations of stochastic logistic growth model in eq. (2) below. The colors blue, orange and green correspond to K=50, 30, 10. In all cases r=1.

To make this model more realistic, let’s see how we can extend it to include stochastic fluctuations (semi-transparent trajectories in the figure above). I’ll look in the following at a simple stochastic logistic growth model (motivated by some discussions with a friend), where the steady state can be calculated exactly. The effect of stochasticity on the steady state is twofold:

  • The population size for long times is not fixed, but fluctuates on a scale \sigma \approx \sqrt{K} around the carrying capacity.
  • The average population size \left\langle n \right\rangle is increased above the carrying capacity K, but the shift goes to 0 as K increases (i.e. the deterministic model is recovered for large K).

Now let’s look at the calculation in detail…

A stochastic, individual-based logistic growth model

In eq. (1), the population is described by a real-valued function N(t). In reality, populations consist of discrete individuals and the population size doesn’t change continuously. So, a more realistic approach is to describe the population dynamics as a birth-death process with the following reactions:

\begin{array}{rl} A & \stackrel{r}{\longrightarrow} A+A \\ A+A & \stackrel{d}{\longrightarrow} A\end{array} (2)

In other words, we assume that during a time interval dt two events may happen:

  • Birth: With a probability r \, dt, any individual may give birth to offspring, thus increasing the population size by 1.
  • Death due to competition: With a probability d\, dt, out of any two individuals one may die due to the competition for common resources. Thus the population size decreases by 1. Note that d is related to the carrying capacity K in eq. (1) by d = r/K.

Note that the stochasticity in this model is not due to random external influences but due to the discreteness of the population (demographic noise).

Solution of the stochastic model

The system of reaction equations (2) translates into the following master equation for the probabilities p_n(t) of the population size at time t being equal to n:

\begin{array}{rl} \displaystyle \partial_t p_n(t) = & \displaystyle n(n+1)d\,p_{n+1}(t) - n(n-1)d\, p_n(t) \\ & \displaystyle + (n-1)r\,p_{n-1}(t) - nr\, p_n(t)\end{array} (3)

This looks daunting. However, it can be simplified a lot by introducing the generating function

\displaystyle G(\lambda,t) := \sum_{n=0}^{\infty} \lambda^n p_n(t)

After some algebra we obtain

\displaystyle \partial_t G(\lambda,t) = \lambda(1-\lambda) \left[d\, \partial_\lambda^2 G(\lambda,t) - r \,\partial_\lambda G(\lambda,t)\right].

This looks simpler than eq. (3) but still finding the full time-dependent solution G(\lambda,t) does not seem feasible. Let’s focus on the steady state where \partial_t G(\lambda,t)=0. Together with the boundary conditions G(\lambda=1,t)=1 and G(\lambda=0,t)=0, we obtain the steady-state solution G_s(\lambda)

\displaystyle G_s(\lambda) = \frac{e^{K\lambda}-1}{e^{K}-1}. (4)

Here we set d=r/K to connect to the notation of eq. (1). Correspondingly, the steady-state probabilities for population size n are

\displaystyle p_n = \frac{K^n}{n!}\frac{1}{e^K-1}. (5)

Of course, these can equivalently also be obtained by solving the recursion relation (3) with \partial_t p_n(t) = 0. This result can be easily checked against simulations, see figure below.

statdist

Stationary probability distribution for the stochastic logistic growth model in eq. (2), for r=1, K=5, d=r/K=0.2. Blue line: exact solution in eq. (5). Yellow bars: Histogram over 10^6 samples.

Steady state with stochasticity

Let’s try to understand what the steady state of our stochastic model looks like. From eq. (4) we easily obtain the first moments of the distribution of the population size. The mean population size is:

\displaystyle \left\langle n \right\rangle = \partial_\lambda\big|_{\lambda=1}G_s(\lambda) = \frac{K}{1-e^{-K}}

So for large carrying capacities K \gg 1, we recover the same result as in the deterministic model (which is reasonable since for large population sizes demographic noise is less pronounced!). However, for smaller K fluctuations increase the average population size. E.g. for K=2, the average population size in our stochastic model is 2/\left(1-e^{-2}\right) \approx 2.31.

Now let us look at the variance,

\begin{array}{rl} \sigma^2 := &\displaystyle \left\langle n^2\right\rangle - \left\langle n \right\rangle^2\\=&\displaystyle\partial_\lambda\big|_{\lambda=1}\lambda \partial_\lambda G_s(\lambda)-\left(\frac{K}{1-e^{-K}}\right)^2\\=& \displaystyle\frac{K\left(1-e^{-K}-K\,e^{-K}\right)}{\left(1-e^{-K}\right)^2}\end{array}.

For large K \gg 1, we have \sigma^2 \approx K. So we see that stochastic fluctuations spread out the steady state to a width \sigma \approx \sqrt{K} around the carrying capacity.

Written by inordinatum

January 1, 2016 at 4:11 pm

How does an absorbing boundary modify the propagator of a stochastic process?

with one comment

A recurring theme throughout my blog are stochastic processes, widely used to model quantities under the influence of both deterministic evolution laws and random fluctuations. A complete solution for a stochastic process is its propagator: The probability density for the final state at some time t_f given an initial condition at time t_i. The propagator depends on the boundary conditions imposed on the stochastic process. For example, an absorbing boundary is useful for describing extinction in a stochastic model of population dynamics: When the population of a given species reaches zero, it cannot be re-introduced.
Typically, however, the “free” propagator for a stochastic process — i.e. the propagator without any boundaries — is much easier to determine than the propagator with any given boundary condition.

However, common sense tells us that the free propagator for a stochastic process contains already all there is to know about it. So I have wondered for quite a while if and how one can recover e.g. the propagator with an absorbing boundary, knowing only the free propagator. I found the solution only quite recently in this preprint on arXiv, and think that it is interesting enough to explain it here briefly.

The result

Consider a one-dimensional, stationary, Markovian stochastic process. Let us denote its free propagator by

\displaystyle G(x_i,x_f|t_f-t_i).

This is the density for the final value x_f at time t_f, given an initial value x_i at time t_i. Stationarity of the process implies that G_0 is a function of the time difference \tau := t_f-t_i only. Define the Laplace transform of the free propagator by

\displaystyle \hat{G}(x_i,x_f|\lambda) := \int_{0}^{\infty} \mathrm{d}\tau\, e^{\lambda \tau}G(x_i,x_f|\tau).

Now let us add an absorbing boundary at x=a. Then, the propagator

\displaystyle G_a(x_i,x_f|t_f-t_i)

is the density for the final value x_f at time t_f, given an initial value x_i at time t_i, without touching the line x=a between t_i and t_f. The upshot of this post is the following simple relationship between the Laplace transforms of the free propagator G and of the propagator with absorbing boundary G_a :

\displaystyle \hat{G}_a(x_i,x_f|\lambda) = \hat{G}(x_i,x_f|\lambda)-\frac{\hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{\hat{G}(a,a|\lambda)}.      (1)

In case the denominator on the right-hand side is zero, one may have to perform some limiting procedure in order to make sense of this formula. We can now check how nicely it works for simple cases, like Brownian motion with drift, where the result is known by other means.

Example: Brownian motion with drift

Let us consider an example, the Brownian motion with drift \mu. Its stochastic differential equation is

\displaystyle  \partial_t X(t) = \mu + \xi(t),

where \xi is white noise with \overline{\xi(t)\xi(t')} =2\sigma\delta(t-t').

The free propagator is a simple Gaussian

\displaystyle  G(x_f,x_i|\tau) = \frac{1}{\sqrt{4\pi\sigma\tau}}\exp\left[-\frac{(x_f-x_i-\mu \tau)^2}{4\sigma \tau}\right].

Its Laplace transform is given by:

\displaystyle  \begin{array}{rl}  \hat{G}(x_f,x_i|\lambda) =& \int_{0}^{\infty} \mathrm{d}\tau\, e^{\lambda \tau}G(x_i,x_f|\tau) \\  =&\frac{1}{\sqrt{\mu^2-4\lambda\sigma}}\exp\left[\frac{\mu}{2\sigma}(x_f-x_i) -\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|x_f-x_i|\right].  \end{array}

To compute G_a, let us assume x_i > a, x_f > a. Note this is without loss of generality, since if x_i and x_f are on opposite sides of the absorbing boundary the propagator G_a is zero. Applying the result (1), we have

\displaystyle   \begin{array}{rl}  \hat{G}_a(x_i,x_f|\lambda) =& \hat{G}(x_i,x_f|\lambda)-\frac{\hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{\hat{G}(a,a|\lambda)}  \\  \quad=& \frac{1}{\sqrt{\mu^2-4\lambda\sigma}}\exp\left[\frac{\mu}{2\sigma}(x_f-x_i) -\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|x_f-x_i|\right] - \\   & \frac{1}{\sqrt{\mu^2-4\lambda\sigma}}\exp\left[\frac{\mu}{2\sigma}(a-x_i) -\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|a-x_i| \right. \\  &\quad \left.+\frac{\mu}{2\sigma}(x_f-a) -\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|x_f-a|\right]  \\  \quad=& \frac{1}{\sqrt{\mu^2-4\lambda\sigma}}e^{\frac{\mu}{2\sigma}(x_f-x_i)}\left[e^{-\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}|x_f-x_i|} - e^{-\frac{1}{2\sigma}\sqrt{\mu^2-4\lambda\sigma}(x_i+x_f-2a)}\right].    \end{array}

The inverse Laplace transform of this expression is:

\displaystyle   \begin{array}{rl}  G_a(x_i,x_f|\tau) =& G(x_i,x_f|\tau)-e^{\frac{\mu}{\sigma}(a-x_i)}G(2a-x_i,x_f|\tau) \\  =& \frac{1}{\sqrt{4\pi\sigma\tau}}e^{-\frac{\mu^2 \tau}{4\sigma} + \mu\frac{x_f-x_i}{2\sigma}}\left[e^{-\frac{(x_f-x_i)^2}{4\sigma\tau}}-e^{-\frac{(2a-x_f+x_i)^2}{4\sigma\tau}}\right].  \end{array}

This is, of course, exactly the same result (the so-called Bachelier-Levy formula) which I had already discussed earlier in a guest blog post at Ricardo’s blog, using slightly different methods (a variant of the Cameron-Martin-Girsanov theorem). Note that some signs which were wrong in the original post are corrected here.

The proof

Now let us go back to the case of a general stochastic process, and discuss how formula (1) can be proven. Let us assume that the stochastic process in question has a path-integral representation, i.e. that

\displaystyle G(x_i,x_f|\tau) = \int_{x(t_i)=x_i}^{x(t_f)=x_f} \mathcal{D}[x(t)]\,e^{-S[x(t)]}

for some action S.

Now, to add an absorbing boundary at x=a we weigh each crossing of x=a by a factor e^{-\gamma} in the path integral, and then let \gamma \to \infty. This effectively retains only those trajectories in the path integral, which never touch x=a. We thus have:

\displaystyle G_a(x_i,x_f|\tau) = \lim_{\gamma \to \infty}\int_{x(t_i)=x_i}^{x(t_f)=x_f} \mathcal{D}[x(t)]\,e^{-S[x(t)]-\gamma\delta(x(t)-a)}.

Now, to evaluate the path integral, we expand it in a series in \gamma (or, in other words, in the number k of crossings of x=a). We get:

\displaystyle   \begin{array}{rl}  \lefteqn{G_a^{(\gamma)}(x_i,x_f|\tau) = G(x_i,x_f|\tau) + }&\\  & + \sum_{k=1}^{\infty} (-\gamma)^k \int_{t_i}^{t_f}\mathrm{d}t_1\, G(x_i,x_1|t_1-t_i) \int_{t_1}^{t_f}\mathrm{d}t_2\, G(x_1,x_2|t_2-t_1) \cdots \\  &  \cdots \int_{t_{k-1}}^{t_f}\mathrm{d}t_k\,G(x_{k-1},x_k|t_k-t_{k-1})\,G(x_k,x_f|t_f-t_k)  \end{array}  .

Note that here I cancelled the factor \frac{1}{k!}, coming from the Taylor expansion of e^{-\gamma \delta(x(t)-a)}, against the k! possibilities of ordering the times t_1\cdots t_k. This is why I can assume them to be in ascending order above. Also, the \delta functions actually set all x_1=x_2=\cdots=x_k=a.
Taking Laplace transforms, the time integrals (which are a k-fold convolution) simplify to a simple product:

\displaystyle  \begin{array}{rl}  \hat{G}_a^{(\gamma)}(x_i,x_f|\lambda) =& \hat{G}(x_i,x_f|\lambda) + \\  & + \sum_{k=1}^{\infty} (-\gamma)^k \hat{G}(x_i,a|\lambda)\hat{G}(a,a|\lambda)^{k-1} \hat{G}(a,x_f|\lambda) \\  =& \hat{G}(x_i,x_f|\lambda) - \frac{\gamma \hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{1+\gamma \hat{G}(a,a|\lambda)}\\  =& \hat{G}(x_i,x_f|\lambda) - \frac{\hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{\gamma^{-1}+ \hat{G}(a,a|\lambda)}  \end{array}.

Now, one easily obtains the claimed result:

\displaystyle  \begin{array}{rl}  \hat{G}_a(x_i,x_f|\lambda) =& \lim_{\gamma \to \infty}\hat{G}_a^{(\gamma)}(x_i,x_f|\lambda) \\  =& \hat{G}(x_i,x_f|\lambda) - \frac{\hat{G}(x_i,a|\lambda)\hat{G}(a,x_f|\lambda)}{\hat{G}(a,a|\lambda)}  \end{array}  .

This computation can, of course, be generalized to more complicated situations like several boundaries, higher dimensions, etc. For details see again the original preprint by C. Grosche. Have fun and let me know if you have questions!

Written by inordinatum

November 24, 2012 at 11:50 pm

A quick introduction to the Martin-Siggia-Rose formalism

with 19 comments

The Martin-Siggia-Rose (MSR) formalism is a method to write stochastic differential equations (SDEs) as a field theory formulated using path integrals. These are familiar from many branches of physics, and are simple to treat perturbatively (or, in some cases, even to solve exactly).

We start from an SDE of the form:

\displaystyle \partial_t x(t) = F(x(t),t) + \xi(x(t),t),

where F is a deterministic function and \xi(x(t),t) is Gaussian noise with an arbitrary correlation function G:

\displaystyle \left\langle \xi(x,t) \xi(x',t') \right\rangle = G(x,t,x',t').

MSR recognized that observables, averaged over solutions of this SDE, can be written as a path integral with a certain action:

\left\langle \mathcal{O}[x(t)] \right\rangle = \int \mathcal{D}[x,\tilde{x}]\,\mathcal{O}[x(t)]\,e^{-S[x,\tilde{x}]},
where
\displaystyle S[x,\tilde{x}] = \int_t i\tilde{x}(t)\left[\partial_t x(t) - F(x(t),t)\right] + \frac{1}{2}\int_{t,t'} G(x(t),t,x(t'),t')\tilde{x}(t)\tilde{x}(t').

A useful special case is an automonous (time-independent) SDE with multiplicative noise:

\displaystyle \partial_t x(t) = F(x) + H(x)\xi(t),

where \overline{\xi(t)\xi(t')} = 2\sigma \delta(t-t') is Gaussian white noise. The resulting action for the path integral is local in time:

\displaystyle S[x,\tilde{x}] = \int_t i\tilde{x}(t)\left[\partial_t x(t) - F(x(t))\right] + \frac{\sigma}{2} H(x(t))^2\tilde{x}(t)^2.

I will now show a few concrete examples, and then explain how this mapping is derived.

Example: Ornstein-Uhlenbeck process

The Ornstein-Uhlenbeck process is one of the few stochastic processes treatable analytically. Its SDE is:
\displaystyle \partial_t x(t) = -\alpha x(t) + \xi(t),
where \alpha > 0 (else the process blows up) and \xi(t) is white noise:
\left\langle \xi(t)\xi(t') \right\rangle = 2\sigma \delta(t-t').

In this section, I will show that the generating functional for this Ornstein-Uhlenbeck process is given by:
\displaystyle \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \exp \left[ \frac{\sigma}{2\alpha}\int_{-\infty}^\infty \mathrm{d} s_1 \int_{-\infty}^\infty \mathrm{d} s_2 \,e^{-\alpha|s_1-s_2|}\, \lambda(s_1) \lambda(s_2) \right] .

To show this, we apply the MSR formula to express the generating functional as a path integral:

\displaystyle \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \int \mathcal{D}[x,\tilde{x}] \exp\left\{-\int_t i\tilde{x}(t)\left[\partial_t x(t) +\alpha x(t)\right] - \sigma\int_{t} \tilde{x}(t)^2 + \int_t \lambda(t) x(t)\right\}.

This path integral can be solved exactly. The key observation is that the exponent is linear in x(t), and hence one can integrate over this field. Using \int_x e^{i x q} \propto \delta(q), this gives (up to a normalization factor \mathcal{N}) a \delta-functional:

\displaystyle \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \mathcal{N}\int \mathcal{D}[\tilde{x}] e^{ \sigma\int_{t} \tilde{x}(t)^2 } \delta\left[i\left(\partial_t - \alpha \right) \tilde{x}(t)+\lambda(t)\right].
The e^{ \sigma\int_{t} \tilde{x}(t)^2 } factor remains since it is the only term in the action independent of x(t).

The \delta-functional fixes \tilde{x}(t) in terms of \lambda:
\displaystyle \tilde{x}(t) = i \int_{t}^\infty \mathrm{d} s \,e^{\alpha(t-s)}\, \lambda(s) .

From this we get
\displaystyle \int_t \tilde{x}(t)^2 = - \int_t \int_{t}^\infty \mathrm{d} s_1 \,e^{\alpha(t-s_1)}\, \lambda(s_1) \int_{t}^\infty \mathrm{d} s_2 \,e^{\alpha(t-s_2)}\, \lambda(s_2) \\  =  - \frac{1}{2\alpha}\int_{-\infty}^\infty \mathrm{d} s_1 \int_{-\infty}^\infty \mathrm{d} s_2 \,e^{-\alpha|s_1-s_2|}\, \lambda(s_1) \lambda(s_2)

The normalization factor \mathcal{N} is seen to be 1 by imposing that for \lambda=0 we must have \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = 1.
Plugging this in gives the generating functional, as claimed above:
\displaystyle \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \exp\left[ \frac{\sigma}{2\alpha}\int_{-\infty}^\infty \mathrm{d} s_1 \int_{-\infty}^\infty \mathrm{d} s_2 \,e^{-\alpha|s_1-s_2|}\, \lambda(s_1) \lambda(s_2) \right] .

One simple consequence is the exponential decay of the connected correlation function:
\displaystyle \left\langle x(t)x(t')\right\rangle^c = \frac{\delta^2}{\delta \lambda(t) \delta \lambda(t')}\big|_{\lambda=0} \ln \left\langle e^{\int_t \lambda(t) x(t)} \right\rangle = \frac{\sigma}{2\alpha}e^{-\alpha|s_1-s_2|}.

But of course, the full generating functional tells us much more – in some sense, it is the exact solution of the SDE!

Note added 01/12/2012: Of course, the generating functional can be written down immediately once one knows the two-point correlation function and the fact that the Ornstein-Uhlenbeck process is Gaussian. However, the latter fact is not immediately obvious (and was not used in our derivation), in fact one can see our MSR calculation as one way to prove Gaussianity for the OU process.

Derivation of the Martin-Siggia-Rose formula

Now that we’ve seen what the MSR path integral formulation is useful for, let us see how it can be derived. The average over an observable can be written as a path integral with the SDE enforced through a \delta-functional:
\displaystyle \left\langle \mathcal{O}[x(t)] \right\rangle = \left\langle\int \mathcal{D}[x]\,\mathcal{O}[x(t)]\,\delta\left[\partial_t x(t) - F(x(t),t) - \xi(x(t),t)\right]\right\rangle.
Now, we rewrite the \delta-functional using the formula \delta(q) \propto \int_{\tilde{x}} e^{i \tilde{x} q}. Since the \delta-functional is a product of a \delta-function at all points in space, effectively we are introducing an auxiliary field \tilde{x}(t):
\displaystyle \left\langle \mathcal{O}[x(t)] \right\rangle = \left\langle\int \mathcal{D}[x,\tilde{x}]\,\mathcal{O}[x(t)]\,\exp\left\{-\int_t i\tilde{x}(t)\left[\partial_t x(t) - F(x(t),t) - \xi(x(t),t)\right]\right\}\right\rangle.
Now, the only thing that depends on the noise \xi(x(t),t) is the factor e^{\int_t i\tilde{x}(t)\xi(x(t),t)}. Since \xi is Gaussian, the average of this factor can be simply expressed as
\displaystyle \left\langle e^{\int_t i\tilde{x}(t)\xi(x(t),t)} \right\rangle = \exp\left[-\frac{1}{2}\int_{t,t'}\tilde{x}(t)\tilde{x}(t')\left\langle\xi(x(t),t)\xi(x(t'),t')\right\rangle\right] =  \exp\left[-\int_{t,t'}\tilde{x}(t)\tilde{x}(t')G(x(t),t,x(t'),t') \right]
Altogether we obtain
\displaystyle \left\langle \mathcal{O}[x(t)] \right\rangle = \int \mathcal{D}[x,\tilde{x}]\,\mathcal{O}[x(t)]\,\exp\left\{-\int_t i\tilde{x}(t)\left[\partial_t x(t) - F(x(t),t)\right]-\int_{t,t'}\tilde{x}(t)\tilde{x}(t')G(x(t),t,x(t'),t') \right\}
which is just the MSR formula claimed above.

Amazingly, I did not yet find a clear derivation of this path integral representation on Wikipedia or anywhere else on the web… So I hope this will be useful for someone starting to work with the formalism!

If there is interest, I may write another blog post in the future with more examples and explanations e.g. on the physical role of the auxiliary field \tilde{x}.

Written by inordinatum

September 27, 2012 at 11:03 pm