Bug in handling events correctly with model changes

[matthiaskoenig]

Hi all,
I just noticed the following bug in roadrunner which results in incorrect model results.
If one has a model with an event trigger and parameters are changed in the trigger before simulation the events do not trigger correctly! I.e. it seems that the event triggers are not correctly evaluated, i.e. not using the actual model values, but the values at model loading.

In the example attached two events are triggered if

IVDOSE_icg > 0 mg

I.e. everytime a simulation starts with IVDOSE_icg > 0 the events should be triggered. I set a value in the SBML model of IVDOSE_icg = 32.5 mg, so the event is triggered if just loading the model.

import roadrunner
import pandas as pd
from matplotlib import pyplot as plt
r: roadrunner.RoadRunner = roadrunner.RoadRunner("icg_body_flat.xml")
r.timeCourseSelections = ["time", "[Cve_icg]", "IVDOSE_total_icg", "IVDOSE_icg"]
s = r.simulate(start=0, end=20, steps=2000)
df = pd.DataFrame(s, columns=s.colnames)
plt.plot(df.time, df["[Cve_icg]"])
plt.show()

And the following curve results.

When setting the value to 0 and back to 32.5 before simulation the model shows incorrect zero!.

import roadrunner
import pandas as pd
from matplotlib import pyplot as plt
r: roadrunner.RoadRunner = roadrunner.RoadRunner("icg_body_flat.xml")
r.timeCourseSelections = ["time", "[Cve_icg]", "IVDOSE_total_icg", "IVDOSE_icg"]

# just setting the value and back ????!!!! This should not change anything, but events are not working any more
r.setValue("init(IVDOSE_icg)", 0.0)
r.setValue("init(IVDOSE_icg)", 32.5)

s = r.simulate(start=0, end=20, steps=2000)
df = pd.DataFrame(s, columns=s.colnames)
plt.plot(df.time, df["[Cve_icg]"])
plt.show()

I have the feeling the underlying problem is that model changes and init changes via setValue do not update the trigger expressions or not reset the trigger expressions (at time zero). I.e. the triggers must be evaluated with the actual values and if a trigger has initalValue=False this must be reset at the beginning of a simulation.

Also when loading the model with IVDOSE_icg = 0 in the SBML file setting the dose via r.setValue("IVDOSE_icg", 32.5) before the simulation does not work to trigger the event. This is how I discovered this originally.

In addition the above example should also work to trigger the event if one sets values in the model, not just the initial values

r.setValue("init(IVDOSE_icg)", 0.0)
r.setValue("init(IVDOSE_icg)", 32.5)

Hope this information helps. But this seems to be a major bug with events in roadrunner.

icg_body_flat.zip

[luciansmith]

The semantics of events are defined in SBML in terms of a simulation, not in terms of user direct manipulation. As I understand it, however, I would assume that events fire when their triggers transition from false to true as part of a simulation. Transitioning from false to true as part of the user manipulating the model state is not really covered, and would require a lot more infrastructure to support. If that would be useful for you, we could put that in as a feature request, but it would, I believe, take a reasonable amount of work. In general, roadrunner watches for trigger transitions during a simulation, not between simulations.

The only exception is for when the model is at its initial state. In that case only, and only when the trigger has a computed value of "true" and an initialValue attribute of "false", will the event be triggered from the beginning of the simulation. In all other situations, it's the computed value that must transition from false to true during simulation.

One thing to consider is the 'oneStep' simulation: one wouldn't want the event to trigger at every single oneStep simulation; only at the first one.

This means that:

r.setValue("init(IVDOSE_icg)", 32.5)

should indeed work as you expect it to. (Setting an 'init' value resets the model. This is a slightly odd design decision, but I didn't want to change it as it would mean a change in behavior if I took it out.) In my greatly simplified example, roadrunner works as I would expect it to in the following program:

import tellurium as te

r = te.loada("""
at S1>1, t0=false: k1 = 1
-> S2; k1
k1 = 0
species S1 = 0
species S2 = 0

""")

r.simulate()
r.plot()
r.setValue("S1", 2)
r.simulate()
r.plot()
r.setValue("init(S1)", 2)
r.simulate()
r.plot()

Does this match the behavior you see? The only time an event can fire 'in an instant' is from the initial state of the model.

[matthiaskoenig]

Yes, the behavior I see is as follows:

import roadrunner
import tellurium as te


r: roadrunner.RoadRunner = te.loada("""
at S1>1, t0=false: k1 = 1
-> S2; k1
k1 = 0
species S1 = 0
species S2 = 0

""")
r.timeCourseSelections = r.timeCourseSelections + ["k1"]

# event not triggered, i.e. k1=0
print("--- c1 ---")
c1 = r.simulate(start=0, end=5, steps=5)
print(c1)

# event not triggered ?, why ?, the state is now S1 = 2; and S1>1 == true => k1 = 1
print("--- c2 ---")
r.setValue("S1", 2)
c2 = r.simulate(start=0, end=5, steps=5)
print(c2)

# event not triggered ?, why ?, the init(state) is now S1 = 2; and S1>1 == true => k1 = 1
print("--- c3 ---")
r.setValue("init(S1)", 2)
c3 = r.simulate(start=0, end=5, steps=5)
print(c3)

results in

--- c1 ---
    time, [S1], [S2], k1
 [[    0,    0,    0,  0],
  [    1,    0,    0,  0],
  [    2,    0,    0,  0],
  [    3,    0,    0,  0],
  [    4,    0,    0,  0],
  [    5,    0,    0,  0]]

--- c2 ---
    time, [S1], [S2], k1
 [[    0,    2,    0,  0],
  [    1,    2,    0,  0],
  [    2,    2,    0,  0],
  [    3,    2,    0,  0],
  [    4,    2,    0,  0],
  [    5,    2,    0,  0]]

--- c3 ---
    time, [S1], [S2], k1
 [[    0,    2,    0,  0],
  [    1,    2,    0,  0],
  [    2,    2,    0,  0],
  [    3,    2,    0,  0],
  [    4,    2,    0,  0],
  [    5,    2,    0,  0]]

I.e. despite the simulation clearly having a trigger change of false (via initValue) to true at the beginning of the simulation at t=0 the event is not triggered correctly and the assignment is not performed k1=1 !? In no way can this be the correct/intended behavior of roadrunner?

The semantics of events are defined in SBML in terms of a simulation, not in terms of user direct manipulation. As I understand it, however, I would assume that events fire when their triggers transition from false to true as part of a simulation. Transitioning from false to true as part of the user manipulating the model state is not really covered, and would require a lot more infrastructure to support.

I don't understand this. Are you saying roadrunner is not handling events correctly because it is too complicated to do this? I.e. the only way to be sure a model with events is executed correctly in roadrunner is to load the model and simulate it, but there is no way to do set a value before simulation? Because the set value affects the state and event handling based on the state is not supported in roadrunner?!
The basic idea of a simulation tool such as roadrunner is to be able to manipulate the model/state and then run simulations with the expected behavior. To have to store every model state with a change as SBML and load it again to be sure it is simulated correctly does not make any sense.

I will test what COPASI is doing in such a situation, but I am sure they handle this as intended, i.e. based on the state. It is even written in the SBML specification that the evaluations are based on the current values of the model!

[matthiaskoenig]

COPASI gives the correct/expected results.
I.e. exporting the model to SBML and using copasi-basico

from basico import *

basico.load_model("tellurium_sim.sbml")
selections = ['Time', '[S1]', '[S2]', 'CN=Root,Model=NoName,Vector=Values[k1],Reference=Value']

print("--- c1 ---")
# tc = run_time_course(0, 5, 5, update_model=True)
tc = run_time_course_with_output(selections, 0, 5, 5)
print(tc)

print("--- c2 ---")
set_species(name='S1', initial_concentration=2)
tc = run_time_course_with_output(selections, 0, 5, 5, update_model=True)
print(tc)

I get

--- c1 ---
   Time  [S1]  [S2]  Values[k1]
0   0.0   0.0   0.0         0.0
1   1.0   0.0   0.0         0.0
2   2.0   0.0   0.0         0.0
3   3.0   0.0   0.0         0.0
4   4.0   0.0   0.0         0.0
5   5.0   0.0   0.0         0.0
--- c2 ---
   Time  [S1]  [S2]  Values[k1]
0   0.0   2.0   0.0         1.0
1   1.0   2.0   1.0         1.0
2   2.0   2.0   2.0         1.0
3   3.0   2.0   3.0         1.0
4   4.0   2.0   4.0         1.0
5   5.0   2.0   5.0         1.0

Changing the init value of S1 so that the event is triggered at t=0 after loading the model triggers the event. COPASI does not allow to set values, but only init values as I understand it so the second case (c3) can not be tested.

*edit (updated COPASI results with full selection)

[matthiaskoenig]

Forgot to comment about the onestep. I see your point here and probably must be handled differently i.e. there should be an option/parameter to decide if the initialValue of triggers is reset for every onestep or not.

I think the main issue is that managing of the init values between simulations with the same model is a bit of a mess (or I have no glue what is happening).

How I expect roadrunner to work (and probably most users) is that at the end of a given simulation the final values are the new init values if a second simulation is run afterwards with the same model. I.e. if there is no reset the second simulation just starts with the values the first simulation left off. I.e. for everything which has an init at the end of the first simulations the init values should be updated in the model with the final values so the next simulation has the correct starting state (and events and other things are evaluated correctly)!. As I understand, this is how other SBML simulators handle this. For instance in copasi-basico the flag exists

`update_model` (bool): sets whether the model should be updated, or reset to initial conditions.

I.e. there are just two options: 1. either reset to initial conditions after the simulation (which are the set of init conditions) or 2. update the model at the end of the simulation so it has the new init conditions the first simulation left off (please correct me if this is not true).

Roadrunner provides a lot of reset functions (see http://sys-bio.github.io/roadrunner/docs-build/PythonAPIReference/cls_ExecutableModel.html?highlight=reset#RoadRunner.ExecutableModel.reset) but none of them is doing the simple update_model, i.e. updating the init conditions of the model at the end of the simulation, so that the next simulation can use the final state as a new init state of the next simulation.

Solution for this issue could be

• to have an option to update_model which at the end of the simulation updates the init values with the final values of the model.
• every call to simulate (and probably also onestep) is a new simulation which starts with the current init values (which would be the final values of the last simulation in case of update_model).
• triggers must be evaluated correctly with the current init values (i.e. if init values were changed by the user or due to update_model this must trigger events, because every simulation is a new simulation with its own t=0)

We could have a quick call. Not sure I make complete sense ;)

[luciansmith]

It looks to me like Copasi doesn't allow you to perform simulations in the mode that roadrunner takes as its default, i.e. 'c2'. It only allows c1 and c3, where it behaves the same as roadrunner. Roadrunner may call the time '0.0', but it has never considered this to be the 'initial state' of the model. To get to the initial state, you need to call 'reset' in some form. Up until the introduction of the 'initialValue' of an event, the differences were minimal to nonexistent!

Again, event triggers simply cannot be triggered in the middle of a simulation without history, because an event is triggered when the model state transitions from false to true. The only exception to this is at the initial state of the model, when an event trigger can be 'true' and its initialState can be 'false'.

But Roadrunner has never considered the state of the model after a simulation to be the 'initial state'. You can simulate for a while, adjust things, then continue the simulation from there. If you adjust things, roadrunner assumes that you are essentially taking over the job of simulator, and any changes that need to be made, you will perform them by hand. Maybe you wanted an event to fire and maybe you didn't; roadrunner has no way of knowing, and leaves those jobs to you.

However, I do see a reasonable use case for what you want. I think I would consider the introduction of a function 'resetCurrentToInit' or something. Essentially, it would take the current state of the model and set all the 'initial' values to those values, and then reset. I could also imagine adding an argument to 'simulate' along the lines of Copasi's 'update_model', where the same thing could happen (set current to initial, set the state of the model to 'initial state'.) Would that work?

[matthiaskoenig]

@luciansmith
So I am getting closer to the problem. Your toy example works as expected. So the issue seems to be with multiple assignments at the same time. I.e. I think the bug is in roadrunner when having dependent assignments in an event. Example file attached.

icg_body_flat.zip

The event has the following assignments

// Events:
  event_IVDOSE_icg_0: at IVDOSE_icg > 0 mg, t0=false: IVDOSE_icg = 0 mg, IVDOSE_total_icg = IVDOSE_total_icg + IVDOSE_icg;

I assume roadrunner is handling the two simultaneous assignments incorrectly. Due to the persistent flag and the simultaneous evaluation of assignments the following assignments should result after setting

init(IVDOSE_icg) = 32

at t=0 in

IVDOSE_icg = 0
IVDOSE_total_icg = 32

which COPASI is doing correctly, but is not working in roadrunner.

import roadrunner
import pandas as pd

r: roadrunner.RoadRunner = roadrunner.RoadRunner("icg_body_flat.xml")
r.timeCourseSelections = ["time", "[Cve_icg]", "IVDOSE_total_icg", "IVDOSE_icg"]

r.setValue("init(IVDOSE_icg)", 32.5)
print(r.getValue("init(IVDOSE_icg)"))

s = r.simulate(start=0, end=10, steps=5)
df = pd.DataFrame(s, columns=s.colnames)
print(df)

results in

32.5
--- roadrunner ---
   time  [Cve_icg]  IVDOSE_total_icg  IVDOSE_icg
0   0.0        0.0               0.0         0.0
1   2.0        0.0               0.0         0.0
2   4.0        0.0               0.0         0.0
3   6.0        0.0               0.0         0.0
4   8.0        0.0               0.0         0.0
5  10.0        0.0               0.0         0.0

Whereas COPASI gives

--- copasi ---
   Time  [Cve_icg]  Values[IVDOSE_total_icg]  Values[IVDOSE_icg]
0   0.0   0.000000              3.250000e+01                 0.0
1   1.0   0.012076              7.948799e-03                 0.0
2   2.0   0.009782              1.944138e-06                 0.0
3   3.0   0.007920              4.762634e-10                 0.0
4   4.0   0.006402              3.546065e-13                 0.0
5   5.0   0.005169             -3.940907e-14                 0.0

created via

# strip SBML names, so COPASI is using ids as names :(
import libsbml
doc: libsbml.SBMLDocument = libsbml.readSBMLFromFile("icg_body_flat.xml")
sbase: libsbml.SBase

for sbase in doc.getListOfAllElements():
    if sbase.isSetName():
        sbase.unsetName()
libsbml.writeSBMLToFile(doc, filename="icg_body_flat_no_names.xml")

from basico import *
basico.load_model("icg_body_flat_no_names.xml")
selections = ['Time', '[Cve_icg]', 'Values[IVDOSE_total_icg]', "Values[IVDOSE_icg]"]

print("--- copasi ---")
# tc = run_time_course(0, 5, 5, update_model=True)

set_parameters(name='IVDOSE_icg', initial_value=32.5)
tc = run_time_course_with_output(selections, 0, 5, 5)
print(tc)

So the bug is in multiple assignments in roadrunner which depend on each other.

[matthiaskoenig]

Further narrowed the issue down and seems only to occur on (multiple) event assignments affecting parameters. Here a minimal example showing the bug. Trying a similar examples with assigning to species seems to work.

import roadrunner
import tellurium as te

r: roadrunner.RoadRunner = te.loada("""
at S1>0, t0=false: S1 = 0, S2 = S2 + S1;

k1 = 1
S1 = 0
S2 = 0

""")
model: roadrunner.ExecutableModel = r.model
sbml = r.getCurrentSBML()
with open("tellurium_sim2.sbml", "w") as f_sbml:
    f_sbml.write(sbml)

r.timeCourseSelections = ["time", "S1", "S2"]

print("--- c1 ---")
c1 = r.simulate(start=0, end=5, steps=5)
print(c1)

print("--- c2 ---")
r.setValue("init(S1)", 2)
c2 = r.simulate(start=0, end=5, steps=5)
print(c2)

--- c1 ---
    time, S1, S2
 [[    0,  0,  0],
  [    1,  0,  0],
  [    2,  0,  0],
  [    3,  0,  0],
  [    4,  0,  0],
  [    5,  0,  0]]

--- c2 ---
    time, S1, S2
 [[    0,  0,  0],
  [    1,  0,  0],
  [    2,  0,  0],
  [    3,  0,  0],
  [    4,  0,  0],
  [    5,  0,  0]]

vs COPASI

--- c1 ---
   Time  Values[S1]  Values[S2]
0   0.0         0.0         0.0
1   1.0         0.0         0.0
2   2.0         0.0         0.0
3   3.0         0.0         0.0
4   4.0         0.0         0.0
5   5.0         0.0         0.0
--- c2 ---
   Time  Values[S1]  Values[S2]
0   0.0         0.0         2.0
1   1.0         0.0         2.0
2   2.0         0.0         2.0
3   3.0         0.0         2.0
4   4.0         0.0         2.0
5   5.0         0.0         2.0

Event is triggered at t=0 and S2 should be 2.0.

[luciansmith]

OK, mucking about with your model, I discovered that the implicit 'reset' that runs when you set an 'init' only resets the species, not everything else. (i.e. it runs 'reset' not 'resetAll'.)

You can fix your example program by changing it to:

import roadrunner
import tellurium as te

r: roadrunner.RoadRunner = te.loada("""
at S1>0, t0=false: S1 = 0, S2 = S2 + S1;

k1 = 1
S1 = 0
S2 = 0

""")
model: roadrunner.ExecutableModel = r.model
sbml = r.getCurrentSBML()
with open("tellurium_sim2.sbml", "w") as f_sbml:
    f_sbml.write(sbml)

r.timeCourseSelections = ["time", "S1", "S2"]

print("--- c1 ---")
c1 = r.simulate(start=0, end=5, steps=5)
print(c1)

print("--- c2 ---")
r.setValue("init(S1)", 2)
r.resetAll() #New line!
c2 = r.simulate(start=0, end=5, steps=5)
print(c2)

There are two obvious fixes on our end: one, to just stop any implicit resets, or two to run an 'resetAll' instead. I'm inclined to go with the latter, as that seems like it matches current behavior a bit better, and then would behave as you expect. Do you have a preference?

[matthiaskoenig]

Thanks for the quick fix. This helps. I don't understand what the implicit resets are doing exactly.
I will test tomorrow with the resetAll if this works for my larger examples.

[luciansmith]

Oh, hey, I figured something out! The thing that tells roadrunner 'I am in the initial state' is whether 'time' has been reset or not. So if you only want to tell roadrunner that it should be in its initial state (without changing any other value), you call:

r.reset(roadrunner.SelectionRecord.TIME)

If you add this to your latest program, you no longer have to call r.setValue("init(S1)'); you only have to set S1 directly:

import roadrunner
import tellurium as te

r: roadrunner.RoadRunner = te.loada("""
at S1>0, t0=false: S1 = 0, S2 = S2 + S1;

k1 = 1
S1 = 0
S2 = 0

""")
model: roadrunner.ExecutableModel = r.model
sbml = r.getCurrentSBML()
with open("tellurium_sim2.sbml", "w") as f_sbml:
    f_sbml.write(sbml)

r.timeCourseSelections = ["time", "S1", "S2"]

print("--- c1 ---")
c1 = r.simulate(start=0, end=5, steps=5)
print(c1)

print("--- c2 ---")
r.setValue("S1", 2)
r.reset(roadrunner.SelectionRecord.TIME)
c2 = r.simulate(start=0, end=5, steps=5)
print(c2)

The whole reason why setting init(S1) worked was solely because of the side-effect of calling 'reset', which resets TIME. So if you instead want to bypass everything else, just reset TIME directly and you're good to go.

[matthiaskoenig]

This sounds great. I hope I understand this correctly.
So instead of setting all current values as init values I could just call the

r.reset(roadrunner.SelectionRecord.TIME)

so that roadrunner uses my current values as init values. I assume this also resets the intialValue of events to the value defined in SBML?! I.e. for roadrunner a second simulation would look like a first simulation if this makes sense.

I.e. I understand this as

# load model
r = RoadRunner(model)
s1 = r.simulate(...)

# => the model has now the values of the final timepoint of the simulation s1

# change values, which overwrite some of the values at the end of s1
r.setValue("S1", 2)

# reset the time => all values in the model are now the new initial values
# initialValue of events has be reset?! to the values when loading the model
r.reset(roadrunner.SelectionRecord.TIME)

# this is a second simulation (which looks like a first simulation after loading the model), but
# with the changes in values from above (simulation s1 and manual changes)
s2 = r.simulate(...)

Is this what happens? If yes, then this is exactly what I am looking for and would solve all my issues. This would avoid the copying of the initial values for the second simulation.

[luciansmith]

That's correct! The only issue is that the values aren't stored as the initial values, so if you later call 'reset' with no arguments, or 'resetAll', it would reset things to the original 'init(x)' values. But if you're not going to do that, you won't need to worry about it.