# ODE solution positivity for systems biology models
**URL:** https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601
**Category:** Modelling & Simulations
**Tags:** differentialequation
**Created:** [March 5, 2025, 11:02pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601 "2025-03-05T23:02:29Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@homegnome0](https://discourse.julialang.org/u/homegnome0)
#### Post date: [March 5, 2025, 11:02pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/1 "2025-03-05T23:02:29Z")
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 6, 2025, 3:53am UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/2 "2025-03-06T03:53:39Z")
Did you try `isoutofdomain` and a low tolerance first?
I haven’t looked into the properties of high hill coefficients in detail but my first thought was they would give positivity.
---
### Author:  [@homegnome0](https://discourse.julialang.org/u/homegnome0)
#### Post date: [March 6, 2025, 5:28am UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/3 "2025-03-06T05:28:08Z")
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.
---
### Author:  [@ranocha](https://discourse.julialang.org/u/ranocha)
#### Post date: [March 6, 2025, 5:46am UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/4 "2025-03-06T05:46:32Z")
Could you please write down the form of the ODE thay causes trouble? It would of course be nice if the model is as simple as possible.
---
### Author:  [@danielwe](https://discourse.julialang.org/u/danielwe)
#### Post date: [March 6, 2025, 6:00am UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/5 "2025-03-06T06:00:24Z")
> [@homegnome0](#):
>
> My thought was that because the RHS is zero when the problem variable becomes zero the solution should always stay positive
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
u(t) = \begin{cases} (t\_\text{Z} - t)^2 / 4 \, , & t \le t\_\text{Z} \, , \\ 0 \, , & t \> t\_\text{Z} \, , \end{cases}
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).
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 6, 2025, 1:56pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/6 "2025-03-06T13:56:51Z")
> [@danielwe](#):
>
> 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
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.
---
### Author:  [@homegnome0](https://discourse.julialang.org/u/homegnome0)
#### Post date: [March 6, 2025, 5:53pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/7 "2025-03-06T17:53:46Z")
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.
---
### Author:  [@homegnome0](https://discourse.julialang.org/u/homegnome0)
#### Post date: [March 6, 2025, 6:00pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/8 "2025-03-06T18:00:22Z")
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.
---
### Author:  [@ranocha](https://discourse.julialang.org/u/ranocha)
#### Post date: [March 6, 2025, 7:09pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/9 "2025-03-06T19:09:46Z")
Thanks! I do not know how the parameters look like, but you may try using some methods from [PositiveIntegrators.jl](https://github.com/SKopecz/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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 6, 2025, 8:25pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/10 "2025-03-06T20:25:43Z")
> [@ranocha](#):
>
> but you may try using some methods from [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl). This package is compatible with [DifferentialEquations.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/DifferentialEquations) and provides some methods guaranteeing non-negative solutions
Totally forgot 😅 yes we should recommend this one more.
---
### Author:  [@homegnome0](https://discourse.julialang.org/u/homegnome0)
#### Post date: [March 6, 2025, 8:55pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/11 "2025-03-06T20:55:38Z")
Thanks, I’ll check it out. Looks like the structure should be just fine since I have exclusively kinetic terms that are positive or negative.
---
### Author:  [@ranocha](https://discourse.julialang.org/u/ranocha)
#### Post date: [March 7, 2025, 5:01am UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/12 "2025-03-07T05:01:55Z")
Please let me know whether it works for you - or whether there are any issues
---
### Author:  [@BdeKoning](https://discourse.julialang.org/u/BdeKoning)
#### Post date: [March 7, 2025, 11:09am UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/13 "2025-03-07T11:09:55Z")
One simple thing you can try is to clip \frac{\text{d}I}{\text{d}t} into [0,\infty) in your RHS when I \< 0.
---
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [March 7, 2025, 11:40am UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/14 "2025-03-07T11:40:44Z")
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:
1. You can just do that coordinate transform
2. 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!)
3. 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:
1. Your solver will care about absolute errors – unless you tell it otherwise via coordinate transforms that transform relative into absolute errors
2. 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.
---
### Author:  [@homegnome0](https://discourse.julialang.org/u/homegnome0)
#### Post date: [March 7, 2025, 6:04pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/15 "2025-03-07T18:04:42Z")
> 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.
---
### Author:  [@homegnome0](https://discourse.julialang.org/u/homegnome0)
#### Post date: [March 7, 2025, 6:05pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/16 "2025-03-07T18:05:45Z")
One issue I forgot: In my full problem I am actually solving DAEs not ODEs, is that supported?
---
### Author:  [@ranocha](https://discourse.julialang.org/u/ranocha)
#### Post date: [March 7, 2025, 6:08pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/17 "2025-03-07T18:08:28Z")
No, PositiveIntegrators.jl does not support DAEs at the moment. This would be an open research topic. How do the DAEs look like?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 7, 2025, 6:31pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/18 "2025-03-07T18:31:15Z")
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)
---
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [March 7, 2025, 6:32pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/19 "2025-03-07T18:32:48Z")
> [@homegnome0](#):
>
> but if the system does reach zero in finite time then these approaches won’t actually help.
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!).
---
### Author:  [@homegnome0](https://discourse.julialang.org/u/homegnome0)
#### Post date: [March 7, 2025, 6:33pm UTC](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601/20 "2025-03-07T18:33:12Z")
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.
[Next page](https://discourse.julialang.org/t/ode-solution-positivity-for-systems-biology-models/126601.md?page=2)