Colored Noise

As discussed in lecture, it may be possible that the noise in a physical or biological system has correlations which are not satisfied by white noise. For example, we found the following equation for the concentration of cGMP:

$$\frac{d [cG]}{dt} = -(\hat{a} n_{PDE}) [cG] + \gamma,$$ (10)

where

$$n_{PDE} = \langle n_{PDE} \rangle + \delta n_{PDE}$$

with

$$\langle \delta n_{PDE} (t) \delta n_{PDE} (s) \rangle = A e^{-k_2 \vert t-s\vert}.$$

This is a stochastic differential equation of the form

$$\frac{dx}{dt} = \alpha x + \gamma + \mu x \xi,$$ (11)

where $\alpha, \gamma$, and $\mu$ are constants, and

$$\langle \xi(t) \xi(s) \rangle = A e^{-\vert t-s\vert/\tau}.$$ (12)

It can be shown that the colored noise $\xi(t)$ may be calculated from the following equation:

$$\frac{d \xi}{dt} = -\frac{1}{\tau} \xi + \frac{\epsilon}{\tau} \eta(t),$$ (13)

where $\eta(t)$ is a Gaussian white noise source with

$$\langle \eta(t) \rangle = 0,$$

$$\langle \eta(t) \eta(s) \rangle = \delta(t-s).$$

Equation (13) describes an Ornstein-Uhlenbeck process, and it can be shown that the solution to this equation satisfies

$$\langle \xi(t) \rangle = 0,$$

$$\langle \xi(t) \xi(s) \rangle = \frac{\epsilon^2}{2 \tau} e^{-\vert t-s\vert/\tau}.$$

Moreover, $\xi(t)$ is Gaussian and stationary only if it is prepared with initial conditions $\xi_0$ consistent with

$$\langle \xi_0 \rangle =0, \qquad \langle \xi_0^2 \rangle = \frac{\epsilon^2}{2 \tau};$$

this means that the initial conditions are chosen from the distribution

$$p(\xi_0) = \left( \frac{\tau}{\epsilon \pi} \right)^{1/2} \exp \left[ -\frac{\xi_0^2}{\epsilon/\tau} \right].$$

Equation (13) can be solved using the Euler-Maruyama method described above, giving $\xi(t)$. This is done with the following program called ou.m:

randn('state',100)
tau = 1; 
xi0 = 1;
dt = 0.01;
N = 1000;
T = N*dt;

dW = sqrt(dt)*randn(1,N);         % Brownian increments
W = cumsum(dW);                   % discretized Brownian path 

xi = zeros(1,N);                 % preallocate for efficiency

xi(1) = xi0 - dt*xi0/tau + dW(1); 

for j=2:N
   xi(j) = xi(j-1) - dt*xi(j-1)/tau + dW(j);
end

plot([0:dt:T],[xi0,xi],'b-'), hold off
xlabel('t','FontSize',12)
ylabel('\xi','FontSize',16,'Rotation',0,'HorizontalAlignment','right')

Text version of this program

This program produces Figure 3.

Figure 3: A solution for the Ornstein-Uhlenbeck process.

Note that in solving for $\xi$ we do not need to know anything about $x$. Equation (11) can then be solved using Euler integration:

$$x_j = x_{j-1} + dt (\alpha x_{j-1} + \gamma + \mu x_{j-1} \xi_{j-1}).$$ (14)

As an exercise, write a program which solves equation (11) using this procedure.