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 ofu(the variables) changing in between. Sometimes one value inubecomes negative somewhere in between the first and the second call tog, the second call toggiving 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