sys-bio / roadrunner

libRoadRunner: A high-performance SBML simulator
http://libroadrunner.org/
Other
39 stars 24 forks source link

Turning on and off conservedMoietyAnalysis doesn't reset total mass. #1146

Open luciansmith opened 11 months ago

luciansmith commented 11 months ago

Here's a short python program that illustrates the problem:

import tellurium as te
r = te.loada("""
      S1 -> S2; k1*S1;
      S2 -> S1; k2*S2;
      S1 = 5; S2 = 10
      k1 = 0.1; k2 = 0.3
""")

print (r.S1, r.S2)
print(r.getSteadyStateValues())

r.conservedMoietyAnalysis = False
r['init(S1)'] = 30
print (r.S2)
print (r.S1, r.S2)
print(r.getSteadyStateValues())
print (r.S1, r.S2)

What's going on in the background is that conservedMoietyAnalysis starts off False by default, and is turned on when we run the model to steady state. We turn it back off explicitly, which lets us set init(S1) without changing S2. But then when we run to steady state again, the conserved moiety analysis turns on again in the background, and it uses the old total mass values ("_CSUM#', I think) instead of re-calculating them.

In addition, the values should be recalculated when setting the current values, too, not just the init:

import tellurium as te
r = te.loada("""
      S1 -> S2; k1*S1;
      S2 -> S1; k2*S2;
      S1 = 5; S2 = 10
      k1 = 0.1; k2 = 0.3
""")

print (r.S1, r.S2)
print(r.getSteadyStateValues())

r.conservedMoietyAnalysis = False
r['S1'] = 30
print (r.S2)
print (r.S1, r.S2)
print(r.getSteadyStateValues())
print (r.S1, r.S2)
adelhpour commented 11 months ago

Could you please try this one:

import tellurium as te
r = te.loada("""
      S1 -> S2; k1*S1;
      S2 -> S1; k2*S2;
      S1 = 5; S2 = 10
      k1 = 0.1; k2 = 0.3
""")

print (r.S1, r.S2)
print(r.getSteadyStateValues())

r._setConservedMoietyAnalysis(False)
r['init(S1)'] = 30
print (r.S2)
print (r.S1, r.S2)
print(r.getSteadyStateValues())
print (r.S1, r.S2)

I am getting different results after setting 'ConservedMoietyAnalysis' to False using this way, but I am not sure if the new results are the ones you are after.

hsauro commented 11 months ago

Don’t compute the steady state values, that might not show up the bug.

luciansmith commented 11 months ago

I get the exact same results with both scripts.

luciansmith commented 11 months ago
5.0 10.0
[11.25  3.75]
10.0
30.0 10.0
[11.25  3.75]
11.25 3.75
luciansmith commented 11 months ago

Note that the problem is that when you change S1, the sum of S1 and S2 is 40, so at steady state, they should total 40. But instead they revert back to total 15.