Simulations with negative start times handle math with csymbol `time` incorrectly

[matthiaskoenig]

Following the discussion today at HARMONY I was trying a few examples to see where I had issues with models with negative start times. One example are events which are triggered based on time>=0. But the problems is basically all math which depends on the csymbol time, which is the relative time to the start of the simulation. The SBML specification states:

An assumption in SBML is that “start time” or “initial time” in a simulation is zero, that is, if t0 is the initial time in the system, t0 = 0. This corresponds to the most common scenario. Initial conditions in SBML take effect at time t = 0. There is no mechanism in SBML for setting the initial time to a value other than 0.

The specification even states that the only way to handle this is to do a time_offset (as mentioned today).

To refer to a different time in a model, one approach is to define a Parameter for a new time variable and use an AssignmentRule in which the assignment expression subtracts a value from the csymbol time.

I.e. the initial time is 0 in a simulation and an event with trigger time>=0 must be triggered at the first simulation point. Same for all triggers and expressions containing the csymbol time; they all refer to the values relative to the start time. This is not the case, but when starting with negative times the event is only triggered if the model time crosses zero.

Example model with an Event:

Event(
    "E1",
    trigger="time>=10",
    assignments={"S_event": 30.0},
)

i.e. at time>=10 the S_event is set. This should happen at 10 after the start time. Another example showing the issue is the model_time which just shows the csymbol time which is by SBML specification 0 at the beginning, but not in the results.

model_t0.zip

import roadrunner
import pandas as pd

r: roadrunner.RoadRunner = roadrunner.RoadRunner("model_t0.xml")
r.timeCourseSelections = r.timeCourseSelections + ["model_time"]
s1 = r.simulate(start=-20, end=20, steps=4)
df1 = pd.DataFrame(s1, columns=s1.colnames)
print(df1)

print()
r: roadrunner.RoadRunner = roadrunner.RoadRunner("model_t0.xml")
r.timeCourseSelections = r.timeCourseSelections + ["model_time"]
s2 = r.simulate(start=0, end=40, steps=4)
df2 = pd.DataFrame(s2, columns=s2.colnames)
print(df2)

Results in:

   time  [S_event]     [S_ia]  model_time
0 -20.0        0.0  20.000000       -20.0
1 -10.0        0.0  18.096692       -10.0
2   0.0        0.0  16.374546         0.0
3  10.0       30.0  14.816335        10.0
4  20.0       30.0  13.406332        20.0

   time  [S_event]     [S_ia]  model_time
0   0.0        0.0  20.000000         0.0
1  10.0       30.0  18.096692        10.0
2  20.0       30.0  16.374512        20.0
3  30.0       30.0  14.816255        30.0
4  40.0       30.0  13.406332        40.0

The [S_event] should be identical for both simulations, as should be model_time.

The correct way of simulation according to the SBML specification is to simulate with start time 0 and time shift the results afterwards. I.e. to simulate

s = r.simulate(start=-20, end=20, steps=4)

the correct way is to do

s = r.simulate(start=0, end=40, steps=4)
s.time = s.time-20

internally. This is basically what is happening in some implementations of SBML solvers, i.e. they just simulate relative to the start time and setting the csymbol time at the beginning to tstart, but this is incorrect. This is why I implemented my simulations all as starting from zero with a timeshift at the end to the desired start time.

COPASI does not even allow to set the start time for a simulation, which is the correct way to do it. Because start is always 0 for SBML models. I.e. the argument: start to r.simulate is incorrect and should probably be renamed to offset or time_shift.

This is also affecting positive start times. I.e. all simulations with start!=0 report incorrect values of the csymbol time.

[fbergmann]

remember that the event also has a flag that can be checked, so that the event will trigger at initial time. Maybe that would take care of it for you.

As for COPASI, there you can specify the models initial time explicitly (though of course that is 0 when importing SBML models). So you could say my current model time is -50, and then simulate for a duration of X. You can then additionally narrow the interval you would like to collect data from, so if you run for a duration of 100, but only want to collect data from 10..50, you could do that.

I think if SBML had a way to specify the models initial time explicitly all those issues would go away.

[matthiaskoenig]

@fbergmann In the COPASI UI I cannot provide a start time. If I add negative values in the "Output Time Points" of the Time Course the negative values are just dropped from the output. I.e. I get the following for -20, -10, 0, 10, 20

Time-Course Task

Problem Description:
    AutomaticStepSize: 0
    StepNumber: 100
    StepSize: 0.01
    Duration: 1
    TimeSeriesRequested: 1
    OutputStartTime: 0
    Output Event: 0
    Start in Steady State: 0
    Use Values: 1
    Values: -20, -10, 0, 10, 20

Method: Deterministic (LSODA)
    Integrate Reduced Model: 0
    Relative Tolerance: 1e-06
    Absolute Tolerance: 1e-12
    Max Internal Steps: 100000
    Max Internal Step Size: 0

Time-Course Result:

# Time  S_ia    Values[model_time]  
0   20  0   
10  18.0967 10  
20  16.3746 20  

I.e. COPASI just drops the negative values

[fbergmann]

@matthiaskoenig here is how to do it in COPASI:

https://share.zight.com/qGuB9dj2