Switching a species between an ode evaluation and a specified function (repeat assignments)

Most of the models in Simbiology will be based upon interacting ODEs. Insulin and glucose, in diabetes models, for example. Insulin affects glucose disposal, and glucose stimulates insulin secretion. In creating, debugging, and fitting such models it is often convenient to "break the loop" and force one of the states to follow a function (of time, for example), to accurately check and calibrate the other system(s). So, one might start by creating the ode model structure in which glucose and insulin are determined by a dynamic mass-balance and rate laws (that is, ODEs). But in some scenarios, it would useful to "switch" the model so that insulin or glucose could be specified as a function over time.

One way to do this would be to create a node that is specified by repeat assignment using a function. In the diabetes example, this would be "functional_glucose". The actual ODE node would be "ode_glucose". One could have a parameter "glucose_switch" that the user would specify (0 for ODE, 1 for functional form). And then, in every rate law that uses glucose, you'd have glucose_switch*functional_glucose + (1-glucose_switch)* ode_glucose. In big models, this would be rather convoluted.

Alternately, you can hook up a repeat assignment rule to the ODE species directly. This avoids extra nodes, and works as long as the "boundary value" box is checked. This is certainly a cleaner way to do this. But if one unchecks the box and runs a simulation, there is an "overspecified model" error thrown because the system is trying to solve the node as an ode, and the repeat assignment rule is fighting the integrator. One could disconnect the repeat assignment rule and uncheck the boundary value box. Is this combination (repeat assignment/boundary value for functional form, and disconnect the assignment rule and uncheck the box for ODE) the best way to have both the ODE and functional form in the model?

In reality, I suspect one would write scripts to hook up the functional form, or to disconnect and let the ODE solver work.

Answer 1

[Arthur Goldsipe]

Hi,

I would take the approach of adding a repeated assignment rule and setting the BoundaryCondition property to true. In some cases, you might just be able to switch between 2 different repeated assignment rules. That's what I did in the following example: Simulating the Glucose-Insulin Response

-Arthur

Answer 2

[Jim Bosley]

Thanks, Arthur.

For others, Arthur's insulindemo example is very helpful: a lot of useful ideas. You can look at how Arthur handled specifying a species value over time.

An example here. Note that if you implement two species with reaction between them, A, B, with stoichiometry A->B with rate law rf = kf*A, and you give A a non-zero initial value, you can easily get the thing to run and get nice exponential traces. If, now, you want to specify A with a function, you have an issue. If you just hook up a repeat assignment rule to A, and (for example) specify A = sin(time), you get an "model overspecified" error. I think you can get rid of the error by changing B to be a boundary node. But then to change back to the ode formulation you have to disconnect the repeat assignment rule, and flip the boundary bit on the node.

What Arthur's insulindemo showed me was that if you set the rate equation stoichiometry to A -> B + A so that no "A" is used Simbiology treats things differently. First, the reaction line from A to the reaction node turns dashed. And now, you can hook up that repeat assignment rule to A and specify A = sin(time) and you get no error. This is interesting, and may be useful, especially if you are doing a network signal model. For mass-balance models I prefer another method.

Let's suppose you have a node inside of a model. Say node "B", where A->B->C. To test your model or to do an experiment you'd like to specify, for example, B = 2 + sin(t). I think it may be better to add a node B_u (for B utilized), and then have two repeat assignment rules setting the value of B_u. One would be active and one not. In one, you specify B_u = B. In the other you set Bu to your desired function, e.g. B_u = 2 + sin(t). In the rate laws between A and B and B and C (and anywhere else B is used in the model) you change the rate law from using B, to use B_u. So, for example, r2 = kf2 * B becomes r2 = kf2 * B_u.

This has some issues, as the mass balance for B can go negative. But if you've removed all mention of B from all rate laws, and are using B_u instead, this should allow model checking, calibration, etc.

I went through Arthur's example and formatted the model so that the graphical representation was available, and made some other minor mods (noted in the .m file) and uploaded it. As I said, the example is excellent and so any errors or inefficiencies introduced are my fault.

Thanks again, Arthur. This was interesting and useful.