stochastic simulations

Equations

The stochastic equations take the following form:

\[\frac{dx}{dt} = P - D + c \sqrt{x} dW \]

where $c$ is the noise strength, $P$ is the production term, $D$ is the degradation term, and $dW$ is the diffusion term modeled as $\mathcal{N}(μ=0,σ=h)$ where h is the step size of the solver.

Numerical integration

We use the Euler-Maruyama solver as implemented in sdeint. The code is produced below for reference.

def noise(x,t):
    c = 4.
    return (c*np.sqrt(abs(x)))

def deltaW(N, m, h):
    # From sdeint implementation
    """Generate sequence of Wiener increments for m independent Wiener
    processes W_j(t) j=0..m-1 for each of N time intervals of length h.    
    Returns:
      dW (array of shape (N, m)): The [n, j] element has the value
      W_j((n+1)*h) - W_j(n*h) 
    """
    return np.random.normal(0.0, h, (N, m))

def eulersde(f,noise,y0,tspan,pars,dW=None):
    N = len(tspan)
    h = (tspan[N-1] - tspan[0])/(N - 1)
    maxtime = tspan[-1]
    # allocate space for result
    d = len(y0)
    y = np.zeros((N+1, d), dtype=type(y0[0]))
    if dW is None:
        # pre-generate Wiener increments (for d independent Wiener processes):
        dW = deltaW(N, d, h)
    y[0] = y0
    currtime = 0
    n = 0
    print(maxtime)
    while currtime < maxtime:
        tn = currtime
        yn = y[n]
        dWn = dW[n,:]
        y[n+1] = yn + f(yn, tn,pars)*h + np.multiply(noise(yn, tn),dWn)
        # Ensure positive terms
        for i in range(len(y[n+1])):
            if y[n+1][i] < 0:
                y[n+1][i] = yn[i]
        currtime += h
        n += 1 
    return y