Re: [SOSlib-discuss] ode solver allows for negative concentrations?

[Rainer M. <ra...@tb...>]

Hi Phillip,

I CCed this to the discuss list, because the question is asked more often.

> I've gotten the sbml ode solver up and running and noticed some curious 
> results. I'm printing out time courses for a simple uni uni reactions 
> with only two species, A and B. By 10 seconds, the output indicates a 
> negative concentrations for species A - which is of course non physical 
> behavior.
> 
> Is there some switch so that I can keep the solver from giving me 
> negative concentrations?

No, there is no such switch.

There are two ways you can get negative concentrations:

1)  The model is wrong:
A correct model of a chemical reaction network model should always obey 
the fundamental law of "mass conservation" and never produce negative 
concentrations!

An irreversible reaction:

A -> B; k * A
dA/dt = -k * A
dB/dt = k * A

If A get's zero, then the rate get's zero and thus A can never decrease 
below 0.

If you choose a kinetic law like e.g. "k", then of course A can get 
negative. However, this doesn't make any sense in chemical kinetics, and 
you would need some justification for such a kinetic law. Such a 
justification could e.g. be the A is assumed to be buffered and always 
constant. In this case A will also never decrease below 0, of course.

Such errors (in mass conservation) are very common in modeling, and thus 
it is actually a frequently used quality check to see if a model 
produces negative concentrations.

2) Numerical instabilities, too high error tolerances
If you are sure that the model is correct, then that's what probably 
happens. You can set very low error tolerances, e.g.

--error 1e-30 --rerror 1e-14

but you should play around with these to find good settings. If the 
solver is getting problems with too low tolerance, then set a higher

--mxstep 1e9

to give the solver more time to solve the problematic region.

Hope that helped,
Rainer