Simulation of SBML models with changing compartment size

[hemilpanchiwala]

This issue is created to discuss the SBML Test Suite models where the values are dependent on the changing compartment size. Currently, the models are failing the tests with deviation in the results. Some of the models (from SBML Test Suite) with this issue are 1504-1514, 1198, 1478, 1498, and some others.

Discussing model 1504, In this model, there is a compartment C1, species S1, and parameter x. The rate rule defined for C1 is 0.2 and for S1 is 0.4 till the event with id _E0 is not triggered.

Now, here as the compartment is changing and the species are in concentration units, the derivative of species S1 will be calculated as below (replacing x by S1 in the image)

And, as the result requires the value of S1 in amount units, so the final S1 value will be

Note: P' -> shows the derivative of P

At t=0, the values are C1 = 0.5, S1 = 0.0, x = 0.0.

At t=1, using the above equation,

So, the value of S1 should be (as the formula is given above)

But the results file from the SBML Test Suite gives S1 = 0.0208 at t = 0.1.

[draeger]

Here is a description how LaTeX code can be rendered within GitHub issues: https://gist.github.com/a-rodin/fef3f543412d6e1ec5b6cf55bf197d7b

An example would be:

[matthiaskoenig]

You cannot calculate the values after t = 0.1. There are additional steps between 0 and 0.1 which have different derivatives and derivatives change between 0 and 0.1 dependent on the intermediate steps. This is what an integration algorithm is doing.

[luciansmith]

There's an sbml-discuss thread from when I created test 1198, which all the other models are variants of:

https://groups.google.com/forum/#!msg/sbml-discuss/jbAFdCkM-fA/l-PANnpmWAgJ

The relevant math bit is:

d(S1_conc)/dt = 0.4 d(C1)/dt = 0.2 S1_conc = S1_amt/C1

You can calculate this a few different ways, the most straightforward of which is simply to calculate S1_conc and C1 from the rates, and then S1_amt from S1_conc and C1. This gives you an exact solution with no uncertainty. If you want to calculate d(S1_amt)/dt from the above equations, you can do that--it gets a little more complicated, but the math works out:

d(S1_amt/C1)/dt = 0.4 (C1dS1_amt/dt - S1_amt dC1/dt) / C1^2 = 0.4 (C1dS1_amt/dt - S1_amt0.2)/C1^2 = 0.4 (dS1_amt/dt) / C1 - 0.2S1_amt/C1^2 = 0.4 (dS1_amt/dt) / C1 = 0.2S1_amt/C1^2 + 0.4 dS1_amt/dt = 0.2S1_amt/C1 + 0.4C1 dS1_amt/dt = 0.2S1_conc + 0.4C1

Feel free to edit this post to make it look more math-y ;-)

[luciansmith]

The upshot is: If the model tells you what the change in concentration over time is, then that is the change in concentration over time: dx/dt is 0.4, and nothing else in the model may change that. (Hence, it changes from 0 to 0.4 over the 1 second of the simulation).

Similarly, dC1/dt is 0.2. This is just a given. In this case, C1 changes from 0.5 to 0.7 over the course of the simulation.

The change in S1's amount over time is a derived quantity based on both C1 and [S1]. When I calculated the canonical answers for this test, I just calculated x*C1 and put the answer in S1's column.

If your simulator has no way to perform these calculations directly, you instead have to figure out dS1_amt/dt in terms of d[S1]/dt and dC1/dt (hence the math). This is a complicated way of figuring out the same thing, but will work out in the end.