In solving a system of biological ODEs, for certain parameter values I am getting solutions that “erroneously” go negative. However, I found this discussion on the DifferentialEquations.jl FAQ page:
“The following set of tools are designed to accuracy enforce positivity in ODE models, which mathematically should be positive in the true solution. If they encounter a model that is actually going negative, they will work really hard to get a positive but correct solution, which is impossible, so they will simply error out. This can be more subtle than you think. Solving u'=-sqrt(u) is not guaranteed to stay positive, even though the derivative goes to zero as u goes to zero (check the analytical solution if you’re curious). Similarly, analyzing nonlinear models can showcase all sorts of behavior. A common cause for accidental negativity is Hill functions in systems biology models: just because derivatives go to zero doesn’t mean they are going to zero fast enough to keep things positive!”
My model indeed uses a Hill function with a coefficient of 4, so I suspect there may be unexpected negative solutions. Can anyone offer some applied math references or examples that could help me determine if this is occurring for my system? Thanks in advance.
I have, which gives me a “Interrupted. Larger maxiters is needed” warning. Increasing maxiters to 1e7 then gives instead “At t=6.0968850103130885, dt was forced below floating point epsilon 8.881784197001252e-16, and step error estimate = 0.0012248955507124634. Aborting. There is either an error in your model specification or the true solution is unstable.” At nominal parameter values, neither of these issues arise and the solutions behave as expected, so the equations themselves shouldn’t have any typos.
My thought was that because the RHS is zero when the problem variable becomes zero the solution should always stay positive, but the FAQ discussion makes me wonder if that’s too naive of an analysis.
Just speaking generally, such a system will remain nonnegative, but not necessarily positive. That is to say, it can reach zero in finite time and with a discontinuous second derivative. Take the u' = -\sqrt{u} example from the docs: the analytical solution is
where the value of t_\text{Z} depends on the initial condition.
Hard to say if this is related to your problem. At least make sure isoutofdomain accepts zero and only rejects strictly negative values. Taking into account truncation error, you should probably expect that small time steps will be needed to stick such a landing without ever going slightly negative, and that comes on top of the small time steps already required to deal with the non-analyticity (ideally you’d set a tstop at any point where a variable hits and sticks to zero, but I assume you don’t know in advance when this may happen).
Shampine has a paper where he argues allowing negativity up to a constant*abstol, because of the natural error at that level effectively being random. It’s a bit difficult to make that work with the implementation of things like sqrt though, so truncated forms could be useful for that.
Certainly, if that would help. The basic form is:
\frac{dT}{dt} = 1 - \theta_{12} T I \\
\frac{dI}{dt} = \theta_{12} T I - E I \\
\frac{dE}{dt} = \theta_3 - \theta_4 E \frac{Q^n}{1+Q^n} \\
\frac{dQ}{dt} = \theta_5 I \\
T(0) = \theta_6 e^{-\theta_{12}} \\
I(0) = \theta_6 \left(1-e^{-\theta_{12}} \right) \\
E(0) = 0 \\
Q(0) = 0.
The variable I is the problem variable that goes negative (or becomes very close to zero) and its log is the main observable.
That is helpful to know and makes more sense then the solution going negative, though it will still cause issues for me because I ultimately take the log of the variable.
isoutofdomain seems to be working correctly, but even before I had it implemented I would find cases where maxiters was exceeded and the integration was interrupted, essentially when the variable was extremely close to zero.
Thanks! I do not know how the parameters look like, but you may try using some methods from PositiveIntegrators.jl. This package is compatible with DifferentialEquations.jl and provides some methods guaranteeing non-negative solutions. You have to adapt your problem definition a bit since these methods require a special structure - but these constraints are satisfied as far as I can tell from the equations you posted above.
I mean, the natural coordinate transform for this kind of thing would be x = log(I), for \dot x = \theta_{12}T-E.
Then you need to decide on your approach:
You can just do that coordinate transform
You can also just add that extra uncoupled variable to your system / integrate x along your trajectories to figure out x(t) (because, as you noticed, integrating I and taking \log I(t) will give messed up answers!)
You can also use adaptive coordinate transforms: Use log coords for small I and non-log coords for large I
The appropriate thing depends on what kind of trajectories you’re looking at, and what kind of observables you care about.
Positivity is just a symptom of your underlying problem, and papering over that via constraints will not help you.
Your underlying problem is that:
Your solver will care about absolute errors – unless you tell it otherwise via coordinate transforms that transform relative into absolute errors
If you care about relative errors, or your dynamics care about relative errors, then your solver’s error estimates / theory won’t mean shit for whatever you care about.
“You care about relative errors” means: log I is the main observable you care about.
“The dynamics care about relative errors” means: Your true trajectories are close to hereroclinics with different stabilities and you care about poincare maps.
The appropriate thing depends on what kind of trajectories you’re looking at, and what kind of observables you care about.
I agree, and this is why my first concern is whether the system itself goes exactly to zero. I realize there many, many different ways of numerically dealing with solutions that come very close, but if the system does reach zero in finite time then these approaches won’t actually help.
You could use MTK to index reduce the DAEs to ODEs, though sending it all the way to index 0 could introduce other numerical errors (specifically in the algebraic constraints, they would no longer be enforced to be exact)
But that’s not your issue! Your right-hand side is smooth / locally lipschitz, Picard-Lindeloef applies, and your convergence rate is at most exponential.
Before your post, there was some discussion about singular systems that reach zero in finite time, like \dot x = - \sqrt{x}. These need special care – and again, positivity is not the main issue, the main issue is that existence of solutions comes from one-sided lipschitz estimates and not from Banach / Picard-Lindeloef, so solvers based on that will obviously suck.
I chalked that up to “unimportant digression” because the system you posted is smooth.
PS: A more instructive example than the square-root thingy is e.g. friction with \dot x \in f(x) with f(x) = -x/|x| for x\neq 0, and f(0)=[-1,1]. All the same considerations apply, but you’re less likely to attempt to use theory that can’t work anyways: You can only solve uniquely in one time-direction, the thing is not smooth, and you need somewhat implicit schemes (which must understand the RHS somewhat, because implicit function theorem / newton iteration won’t solve your nonlinear equation for each time-step!).
For the current case, they are so trivial that they can be written as explicit ODEs, but ultimately the applications will involve general implicit nonlinear functions.
yes we should recommend this one more.