must I use events for discontinuous functions?

A previous question asked about input a function of time. A function y = exp(-k*t) was input, with errors. The answer from Mathworks was that Simbiology was designed for ODEs, and recommended putting the expression in ODE form. I hadn't read that, and put in a rate expression (as a repeat assignment to rate r1, which was used in a rate law) of r1 = alpha * time * exp(-time/beta).

Simbiology may have improved since then, though, because this worked when I used the variable "time". However, I want to repeat this function periodically in simulations. I don't see how dosing could be used for an expression of time. Not sure how events would work, either. I want to repeat this (for example, every 1440 minutes). So I assumed that "time" was a system variable and used it in a repeat assignment of clocktime = mod(time,minperday), with minperday=1440 minutes. This gave results that were not expected. clocktime ended up being delivered as a piecewise continuous function that changed with a piecewise constant rate, with discontinuities somewhat random in timing and magnitude. I couldn't figure out how SB was getting the strange shape.

So I've created a species falsetime. Because it's a species, I had to use units of grams. d(falsetime)/dt = timefont. timefont = 1 gram/minute. Then I used liecon (converting the lie) of 1 minute/gram. And so I can then get clocktime = mod(falsetime*liecon, minperday) to get a sharktooth function. Problem is, the integrator is not good at picking up discontinuities (not criticizing, most extrapolation integrators are bad at this). So the sharktooth function is somewhat irregular. To fix this, I set the max time step in the integrator to 1 minute. This gives me a fairly uniform shark tooth. The question arises, could I set up an arbitrary event for every 1440 minutes, to force the integrator to catch my discontinuity? Is there a better way to do what I'm trying to do?

Answer 1

[Arthur Goldsipe]

Hi Jim,

There are several issues in your post that I'd like to address. I'll try to go in the order they appear:

(1) must I use events for discontinuous functions?

In order to get accurate simulation results, you must use either doses or events to tell the solver when discontinuities occur. I put together this simple example on the topic.

(2) What's wrong with using a mathematical expression like y = exp(-k*t)?

I assume you're referring to this MATLAB Answers post? The answer is incorrect (or at least very misleading), and I will work with Technical Support to update it. The correct answer is that you can use expressions like these in repeated assignment rules, but you must use time instead of t when you want to refer to simulation time. Also note that MATLAB is case sensitive, so you cannot use Time or TIME.

There's a secondary question of when is it appropriate to use repeated assignment rules versus rate rules versus reactions. All are supported with SimBiology, and there are a lot of factors that might affect which you use. But here's my rule of thumb: If I can easily write my model using repeated assignment rules, I do so. If not, I use reactions for species and rate rules for everything else.

I prefer repeated assignment rules because they are faster and more accurate than solving differential equations. I prefer reactions over rate rules because they simplify the bookkeeping and more closely match how I think of the model.

(3) What's going on when I use a repeated assignment rule for clocktime = mod(time,minperday)?

I suspect the plot looks strange because the ODE solver is sampling your "shark tooth function" at irregular times. Any plot of a periodic function will be extremely sensitive to the sampling times. For example, if you plot cos(time) for time=0:2*pi:100*pi it will look like cosine always has a value of 1.

This problem is exacerbated by the fact that your function is both periodic and discontinuous. To get a plot that shows the true shape of the curve, you need to sample just before and just after every discontinuity. Unless the ODE solver knows where those discontinuities occur, it will not be able to sample times appropriately.

(4) What's the best way to model a "shark tooth function" in SimBiology?

This is the crux of your post, but it's important to understand the background above before I provide an answer. My recommendation is to continue implementing the "shark tooth function" using a repeated assignment clocktime = mod(time,minperday). However, you also need to inform the solver about the periodic discontinuities using an event or dose. This is quite similar to the example I mentioned earlier. Since your discontinuities are periodic, I think the easiest approach is to use add a "dummy" repeat dose targeting an arbitrary species in your model. Set the dose amount to 0 so that it doesn't really modify the species, set the interval to 1440 minutes, and the repeat count to a value large enough to cover your entire simulation. The net result is that this dose will force the solver to reset its internal state every 1440 minutes so that it can accurately capture the behavior of your "shark tooth function."

I hope this covers everything.

-Arthur

Answer 2

[Jim Bosley]

Arthur, this is spot on. I had gotten pretty accurate results by setting the max step-size to 1, but that did lower the numerical efficiency in integration a bit (that's meant to be understated humor, btw). The arbitrary zero dose every 1440 is brilliant. I'm simulating a human: that dose can be a tablet of fish oil every day. :)

The post you referred to was the one I'd seen. So using "time" is the key. Are there any other reserved words in simbiology expressions, other than matlab functions? Is there a list somewhere?

Your points on the philosophy (of repeat assignments vs reactions) is very helpful. When there are many ways to go, in a complex situation, an expert's advice to "here's how to start" is invaluable.

The sampling time issue, along with the aliasing example you cite, may be the issue. I'll see if I can put together a clean example that shows the problem (and I'll try to use Shannon's sampling theorem to enforce a meaningful graph. I'll repost any learnings here (that is, I'll confirm a real issue, or I'll confirm that this was a sampling/aliasing issue.

Your cogent response covers everything indeed, many thanks! Jim

Answer 3

[Arthur Goldsipe]

Hi,

The only other reserved word is "null", which is used in reactions when there are no reactants or no products. For example, "null -> x" denotes a production reaction for x. The restrictions on names are discussed here: https://www.mathworks.com/help/simbio/ref/name.html