When solving non-diagonal noise SDE: Function g called twice in each timestep, changing u inbetween?

In short: When making solving SDEs the noise function seems to be called twice in every timestep, with the variable values being changed in between. This causes out of domain errors when my positive_domain callback cannot ensure that all values are positive between the two calls to the noise function. First I describe the code I used when encountering the problem, next a more minimal example.

The origin of this problem is that I try to make simulations of a biochemical reaction network, having a positive_domain callback which in each timestep checks whenever my solutions have become negative and fixes that (in case noise makes the solution negative, which might cause a crash during the noise calculation when the method takes sqrt of something negative).

In principle my method looks something like this:

using DifferentialEquations

kBw = 3600; kDw = 18; kD = 18
kB1 = 3600; kB2 = 3600; kB3 = 3600; kB4 = 1800; kB5 = 3600;
kD1 = 18; kD2 = 18; kD3 = 18; kD4 = 1800; kD5 = 18;
kK1 = 36; kK2 = 36;
kP = 180; kDeg = 0.7;
v0 = 0.4; F = 30; K = 0.2;
λW = 4; λV = 4.5;
η = 0.101
p = [kBw, kDw, kD, kB1, kB2, kB3, kB4, kB5, kD1, kD2, kD3, kD4, kD5, kK1, kK2, kP, kDeg, v0, F, K, λW, λV, η]
@reaction_func vB(σB) = v0*((1+F*σB)/(K+σB))
my_model = @reaction_network rnType η begin
    kDeg,       (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅
    kDeg,       vPp ⟶ phos
    (kBw,kDw),  2w ⟷ w2
    (kB1,kD1),  w2 + v ⟷ w2v
    (kB2,kD2),  w2v + v ⟷ w2v2
    kK1,        w2v ⟶ w2 + vP
    kK2,        w2v2 ⟶ w2v + vP
    (kB3,kD3),  w2 + σB ⟷ w2σB
    (kB4,kD4),  w2σB + v ⟷ w2v + σB
    (kB5,kD5),  vP + phos ⟷ vPp
    kP,         vPp ⟶ v + phos
    vB(σB),     ∅ ⟶ σB
    λW*vB(σB),  ∅ ⟶ w
    λV*vB(σB),  ∅ ⟶ v
end kBw kDw kD kB1 kB2 kB3 kB4 kB5 kD1 kD2 kD3 kD4 kD5 kK1 kK2 kP kDeg v0 F K λW λV

function positive_domain()
    function condition(u,t,integrator)
        minimum(u) < 0
    end
    function affect!(integrator)
        integrator.u .= integrator.uprev
    end
    return DiscreteCallback(condition,affect!,save_positions = (true,true)) 
end
function equi(model, p)
    prob = ODEProblem(model,fill(1.,length(model.syms)),(0.,100.),p);
    sol = solve(prob,Rosenbrock23())
    return max.(sol[end],0.0)
end;
function stochsim(model, p, tspan)
    u0 = equiN(model,p)
    prob_sde = SDEProblem(model,u0,tspan,p)
    sol = solve(prob_sde,ImplicitEM(),callback=positive_domain(),saveat=0.001);
end

p[end] = 0.001
sol = stochsim(my_model,p,(0.,2));

where p[end] = 0.001 is a factor linearly scaling the amount of noise in the model. Here is where the weird stuff is starting. If I run the code it works fine if I set p[end] = 0.001 or p[end] = 0.1. However, if I set p[end] = 0.01 I receive DomainError (not every time, but about every other time I run it). I have traced this error to some variable being pushed into negative value by noise, and then the square root of this value is taken when the amount of noise is calculated. The square root of the negative number yields the error. The problem is that the positive_domain callback should catch this.

Some investigation gives that this is what happens:

  • The callback is called at the beginning of every time step, ensuring that the value is non-negative.
  • In every time step the noise function, g, is called twice, with the value of u (the variables) changing in between. Sometimes one value in u becomes negative somewhere in between the first and the second call to g, the second call to g giving the out of domain error.

Ignoring the case where I found the problem and looking at a minimal example I can use this code:

using DifferentialEquations
function f(du,u,p,t)
    du .= 1.01u
end
function g(du,u,p,t)
  println("Calling g at ",t)
  println(u,"\n")
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))
sol = solve(prob,ImplicitEM());

to confirm that g is called twice at every timestep, with the value of u changed in between. Thereby causing problems is one value in u is made negative after the first g call, but without the positive_domain callback having a chance to fix it before g is called again.

Is this supposed to be like this? If so, is there a way around it to prevent the domain error?

Cheers

File an issue. This’ll get lost here. I’m not sure if it’s an issue or not but won’t be able to look into it in details until a flight tomorrow.

Thank you.

The issue can be found at https://github.com/JuliaDiffEq/StochasticDiffEq.jl/issues/106 in case anyone stumbles upon this thread in the future while looking for answers and want to learn more.