What's the best way of modelling changing volumes?

[julibeg]

Hi there, I just started out with tellurium and it's looking great!

One thing I can't quite figure out though is how to best deal with changing compartment volumes. For example, if I wanted to simulate decreasing concentrations of my species due to dilution because of the compartment volume increasing, how would I best go about that? I know that there is the substanceOnly keyword and that I can define the volume with an assignment rule and time, but was wondering if there is a more canonical way to represent changing volumes. I couldn't really find anything tangible in the documentation. Any hints (or ideally a simple working example) would be greatly appreciated.

Thanks!

[luciansmith]

You can change a compartment volume just like you change any other variable! Assignment rules, rate rules, and events all work.

In this particular case, it sounds like you probably want a rate rule. In Antimony, that's just something like:

C' = 1.3

to indicate that dC/dt is 1.3 (any formula is fine there; it doesn't have to be a constant).

The deal with 'substanceOnly' is that normally, when a species ID ('X') shows up in math, that ID means 'the concentration of X'. However, if you define X as a 'substanceOnly' species, that means that when it shows up in math, that ID means 'the amount of X'. In a model with changing volumes, it's obviously very important to get this right ;-) If, say, your reaction rates increase with an increase of the concentration of X, and your rate law is something like 'k1X', X should be a normal non-substance-only species. If your rate law instead is "k1X/C', then that's when you want X to be the substanceOnly species.

Here's a simple worked example:

import tellurium as te

r = te.loada("""
   S1 -> S2; k1*S1;
   S1 = 10
   S1 in C1
   S2 in C1
   C1 = 1
   C1' = .4
   k1 = 0.4
   S2 = 0

""")

r.simulate()
r.plot()

In the plot, you can see both concentrations slowly go down in the end as the volume continues to increase. If you were to plot S1 and S2 amounts, instead, you wouldn't see the decrease (but would still see the decrease in the rate at which S1 is converted to S2.

[julibeg]

Ah, nice! I didn't realise that one can simply change the volume of a compartment in this way (in my workaround I defined a volume variable and used it to create concentration variables from the substanceOnly species which would then go into the rates). This definitely makes the model definition more straightforward than my original approach. Thanks!

Perhaps adding a short section (e.g. "How to work with changing volumes") to the docs might be helpful for other beginners who are not familiar with Antimony yet.

[matthiaskoenig]

@julibeg Yes, compartment sizes can be changed dynamically.
Just be careful if you use other SBML tools beside tellurium/roadrunner for analysis. Many of them do not handle changing compartment volumes (e.g. AMICI and other tools for parameter fitting, JWS online, ...). If you only use roadrunner/tellurium for simulation (or COPASI) things will work fine.

[julibeg]

You can change a compartment volume just like you change any other variable! Assignment rules, rate rules, and events all work.

@luciansmith sorry for being dense, but would you mind showing how you would use an Assignment Rule on a compartment? I tried to modify your snippet from above, but the below (and other variations) all didn't work

r = te.loada(
    """
    S1 -> S2; k1*S1;
    S1 = 10
    S1 in C1
    S2 in C1
    C1 := 1 + 0.4 * time
    k1 = 0.4
    """
)

print(r.simulate(0, 5, 6, selections=["time", "C1", "S1", "S2"]))
# this prints
#     time,  C1, S1, S2
#  [[    0,   1, -0, -0],
#   [    1, 1.4,  0,  0],
#   [    2, 1.8,  0,  0],
#   [    3, 2.2,  0,  0],
#   [    4, 2.6,  0,  0],
#   [    5,   3,  0,  0]]
# C1 increases, but S1 and S2 are 0

[luciansmith]

That's a bug! Somehow tellurium is getting confused and thinking that S1 is initialized to 0. You can work around the bug with by setting S1 back to 10 expliciitly, i.e.:

import tellurium as te
r = te.loada(
    """
    S1 -> S2; k1*S1;
    S1 = 10
    S1 in C1
    S2 in C1
    C1 := 1 + 0.4 * time
    k1 = 0.4
    """
)
r.S1 = 10
print(r.simulate(0, 5, 6, selections=["time", "C1", "S1", "S2"]))

My guess is that something about putting 'time' into an assignment rule is causing an unexpected result. I'll look into it.