PositiveDomain on all ODEs/SDEs?

[ChrisRackauckas]

@isaacsas had an interesting point that since concentrations and amounts should always be positive, why not add PositiveDomain on them by default? @TorkelE you had some example problems that were hitting 0: what happens if we try that on your problems?

https://github.com/JuliaDiffEq/DiffEqBiological.jl/pull/87#issuecomment-460848662

[TorkelE]

So I am starting to look at this, but have been updating to julia 1.1.0 and the latest version of DifferentialEquations so. This have caused some issues.

Lets start with some code. Importing the required things:

using DifferentialEquations
using DiffEqBiological
using Plots

the model:

#Setting parameter values.
v0 = 0.005; v = 0.1; K = 2.8; n = 4;
kD = 10; kB = 100; k = 0.1;
deg = 0.01; i = 0.; 
η = 0.3;

#The model and the parameter set
pars = [v0, v, K, n, kD, kB, k, deg, i, η]; 
model = @reaction_network rnType η begin
    v0 + hill(X,v,K,n), ∅ → (X+Y)
    deg, (X,Y,XY) → ∅
    (kB,kD), X + Y ↔ XY
    i*k, XY → X
end v0 v K n kD kB k deg i;

and then some utility functions:

#Finds a steady state (ish) to start the simulations from.
function equi(model,p)
    uInit = fill(1., length(model.syms))
    sol = solve(ODEProblem(model,uInit,(0.,100.),p),Rosenbrock23());
    return max.(sol[end],0.0)
end
#Make a single SDE simulations from the model
function stochsim(tspan, model, p;  step_t=-1., step_amp=0., saveat=0.1,adaptive=true,dt=0.001)
    u0 = equi(model,p)
    prob = SDEProblem(model,u0,tspan,deepcopy(p))
    return solve(prob,ImplicitEM(),callback=step_cb(step_t,step_amp),saveat=saveat,adaptive=adaptive,dt=dt,tstops=[step_t])
end
#At a given time, increase the input parameter to a desired value.
function step_cb(t,amp)
    condition(u,t,integrator) = (t==t)
    affect!(integrator) = integrator.p[9] += amp
    return DiscreteCallback(condition,affect!,save_positions = (true,true)) 
end;

(The system will typically start of in an inactive state. This will feature low concentrations routinely dipping into negatives. At a given time we add the input, which will activate the system and move it into larger concentrations. In this case the time from input to actual activation is highly stochastic and may vary widely between different simulations)

Now I initially just wanted to test running this without some kind of positive domain callback, just testing with abs(...) added to the noise term. I now run the SDE:

sol = stochsim((0.,1000.), model, pars)

this yields the following error:

┌ Warning: dt <= dtmin. Aborting. There is either an your model specification or the true solution is unstable.
└ @ DiffEqBase /home/Torkel.Loman/.julia/packages/DiffEqBase/PvfXM/src/integrator_interface.jl:156

looking at the solution using plot(sol) can confirm instability as it suddenly blows up. I have tried for the various stiff SDE solvers.

Next, to try and see what effect the abs() term had I remove that from DiffEqBiological (making line 344 in reaction_network.jl ...sqrt($(deepcopy(reactions[i].rate_DE))))). To prevent going into negatives I make this simple modification

function stochsim_pd(tspan, model, p; step_t=-1., step_amp=0., saveat=0.1,adaptive=true,dt=0.001)
    u0 = equi(model,p)
    prob = SDEProblem(model,u0,tspan,deepcopy(p))
    return solve(prob,ImplicitEM(),callback=CallbackSet(positive_domain(),step_cb(step_t,step_amp)),saveat=saveat,adaptive=adaptive,dt=dt,tstops=[step_t])
end
function positive_domain()
    condition(u,t,integrator) = minimum(u) < 0
    affect!(integrator) = integrator.u .= integrator.uprev
    return DiscreteCallback(condition,affect!,save_positions = (false,false)) 
end; 

(Admittedly not the proper positive domain callback. I've forgotten why I started using this one, maybe it was related to performance and this one working good enough. Or just not getting the main one to work properly. Was about to do some comparison with the various approaches when trouble started happening)

Rerunning sol = stochsim((0.,1000.), model, pars) gives the following error

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase /home/Torkel.Loman/.julia/packages/DiffEqBase/PvfXM/src/integrator_interface.jl:150

which is a different one compared to the previous. There is definitely something fishy going on, the second fault suggests that at least parts of this is me not handing the updated DifferentialEquations properly.

More general I think it is sensible to ensure that the SDE does not go into negative. Unless there's performance implications making PositiveDomain default might be the most sensible option. However, I remember there are some cases where the adaptive time step estimates step sizes, and thus calls the noise function between steps. This cause a call to square root of negative despite the presence of a positive domain callback, so it might be sensible to have the abs(..) there in any case. I do not think it should have a significant implications on performance (or anything else).

[isaacsas]

I noticed the maxiters warning is also showing up in DiffEqBiological tests now, from "func_test.jl", when solve is called on each of the three SDEProblems. I'm not sure what it is referring to since the test uses fixed timestep EM.

[isaacsas]

On how to handle solutions going negative; my understanding is that the underlying CLE models are not really well-defined once a solution component becomes zero. So I'm not sure one should trust any results obtained after that time. Maybe we should provide a callback to back out this time and then terminate the simulation?

I guess it could be made robust in a numerical sense with the abs and positive domain callback, but is that then just misleading users into trusting simulation results whose meaning is dubious?

[TorkelE]

I guess it is up to how much trust one have in that the user keeps track on how things work.

E.g. in this case (or at least in some similar one) I am mostly interested in what happens after an input (where the approximation is better), but it is convenient to start of the run in low concentrations (where the approximation is less correct, but I care less about this incorrectness).

[ChrisRackauckas]

My favorite demonstration of what's going on is u' = -sqrt(u), u(0) = 1. That seems like a simple equation, but actually it has no solution after a finite time.

https://www.wolframalpha.com/input/?i=u%27+%3D+-sqrt(u),+u(0)+%3D+1

The solution only exists for a finite amount of time when you have square roots because it approaches the 0 non-asymptotically: it doesn't slow down fast enough. Here you have the SDE capturing a feature of the real discrete stochastic equation. When you have small copy numbers, you won't go asymptotically to zero: you will hit zero. So these models hit zero, and by numerical error are slightly above or below it. Above it is fine. Below it causes problems.

And it's not even just numerical that it hits zero: the true strong solution in many of these cases will actually hit zero. When that happens, without extra information, there no longer exists a strong solution. So IIRC du = -u dt + sqrt(u)dW_t is an SDE without a strong solution (or something like it). But we can make there be a solution by defining which possible branch we want to take at zero (this would need to be proved, but go with it). We can just throw a physical constraint on there that <0 means 0 with numerical error. What it's kind of saying is that our copy numbers are low enough that the SDE is not a great approximation to the discrete equation, but we can ignore that and still use it under this approximation.

So I think that we should explicitly zero-clamp everything. The easiest way to do this is to always use max(u,(0,)) and then put a callback in that, if <0, -> 0. That technically introduces discontinuities at zero, so the analysis would be difficult, but in some sense it should help with model breakdown. A more intense fix would be to transform all of the variables. Instead of using u, we can write the differential equation on u^2 or exp(u) using Ito's rule. Then we can take in the user's array, transform it, and then run the ODE in that space. The solutions that they will be given are the transformed variables, so to plot correctly they'd want to transform back. But we can add some idea of domain transforms to the SDEProblem via its trait and plug that into the other stuff. Or let the user wrap it themselves. This is of course a lot more difficult but would be a good way to handle it.

[isaacsas]

Another fix, which would also be a fair amount of work, would be to provide a hybrid algorithm that dynamically partitions reactions involving species with small populations to be represented as jumps. I know there are is at least one method for partitioning / coupling as

jump <-> tau-leap <-> CLE

The bigger issue is that there are lots of moving parts to making something like that work...

[ChrisRackauckas]

That's what it's all building towards though. That's why I made sure the jumps work in ODEs and SDEs. And I think going straight to ODEs might be better if you already have tau-leap. Really you just need a mutable array of UInt8's that say where that value is being calculated at right now, and then handle it as a jump or not based on that.

[isaacsas]

Sure, I'm in complete agreement that is something we should build towards. I'd like to be able to use that functionality myself (actually, to use it in spatial systems...).

At least on the jump end, we'll need good functionality for dynamically adding/removing jumps from the SSAs (or perhaps easier, to add support for ignoring them through the mutable array you mention).