Events which introduce dicontinuities (Simbiology toolbox)

Hi,

What is the best way of using events in SimBiology which change the parameters in a discontinuous way. Sometimes the solver crashes because of these which is annoying. I want to simulate an experiment in which the medium is changed from glucose -> glycerol , and back, many times. So the concentrations of these are not continuous, which gives me troubles sometimes. I use sundials solver.

Best,

Razvan

Answer 1

[Arthur Goldsipe]

Hi Razvan,

There's a SimBiology demo on this topic. Execute the following from the MATLAB prompt and see if that answers your question:

showdemo('discontSimBiologyModel')

If you still have questions after reading the demo, please post again.

Also, can you clarify what you mean by "the solver crashes"? If you introduce a discontinuity in an ODE without using events, I would not be surprised if you got some sort of error message that the solver that the step size could not be reduced enough to meet the error tolerances. But if you got some other sort of error, then I encourage you to give further details (reproduction steps, etc.), either to technical support or via email. (I will separately send you an email address for such feedback.)

The reason you may receive an error message is that the solver tries to use a continuous interpolating function to predict the system's behavior across the discontinuity. It sees that the interpolant works poorly and keeps trying to improve the interpolant by reducing the size of the time step that it takes. If the solver reaches the minimum step size and is still unable to ensure the solution's accuracy, it will error. By using an event, you explicitly inform the solver there is a discontinuity, and it uses different interpolants before and after the event.

-Arthur

[Razvan]

I want to mention that Arthur also indicated in an email that I should check the Michalis Menten rates which create instabilities. When the substrate becomes close to zero, because of the absolute tolerance used by the ODE solver, sometimes some concentrations appear to be negative. These can become even more negative. In order to stabilize the equations with respect to negative concentrations one workaround is to change the reaction rates like this:

Instead of having a rate V * [E] * [S] / ( K + [S] )

Use V * [E] * max(0,[S]) / ( K + [S] )

which will give a zero rate for negative [S] (which actually means [S]=0 +/- some error).

If this is not enough, I guess you can use this trick - min(0,[C]) instead of [C] - for all concentrations which you believe to go close to zero.

Thanks again Arthur!