Parameter changes which affect compartment volumes must update species concentrations/amounts

[matthiaskoenig]

Parameter changes which affect compartment volumes must update the concentrations of the species in the compartments. This is currently not done in AMICI resulting in incorrect initial concentrations/amounts and incorrect forward simulations.

Example attached:

"""AMICI example simulation"""
from pathlib import Path
import amici
import numpy as np
import pandas as pd

mid = "icg_sd"
base_path: Path = Path(__file__).parent
model_dir = base_path / mid
model_path = base_path / f"{mid}.xml"

sbml_importer = amici.SbmlImporter(model_path)
sbml_importer.sbml2amici(mid, model_dir)

# load the model module
model_module = amici.import_model_module(mid, model_dir)
model = model_module.getModel()

# instantiate solver
solver = model.getSolver()

# update BW parameter which changes volumes via rules!
model.setParameterById("BW", 80)

# set timepoints
timepoints = np.linspace(0, 10, num=11)
model.setTimepoints(timepoints)

# simulation
rdata = amici.runAmiciSimulation(model, solver)
xids = model.getStateIds()
df = pd.DataFrame(rdata.x, columns=xids)
df.insert(loc=0, column="time", value=timepoints)
print(df)

This results in the following concentrations:

    time  Cre_plasma_icg  Cgi_plasma_icg  Cli_plasma_icg  Clu_plasma_icg  Cve_icg  Car_icg  Cpo_icg  Chv_icg  Afeces_icg  LI__bil_ext  LI__icg  LI__icg_bi  IVDOSE_icg  cum_dose_icg
0    0.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
1    1.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
2    2.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
3    3.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
4    4.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
5    5.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
6    6.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
7    7.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
8    8.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
9    9.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0
10  10.0             0.0             0.0             0.0             0.0      0.0      0.0      0.0      0.0         0.0         0.01      0.0         0.0         0.0           0.0

Because the compartment volume of LI__bil_ext was changed the concentration has to change. For the given BW change from 0.01 to 0.094.

See comparison of the simulation with roadrunner, COPASI gives the identical results as roadrunner, i.e. the 0.009375. By changing the BW from 75 to 80 the volume changes by 80/75 and the concentration by 0.01 * 75/80 = 0.009375

Model attached: icg_sd.zip

Probably the same bug happens when updating compartment volumes. This is a major bug and all changes resulting in compartment volume changes will have incorrect initial states.

[FFroehlich]

This definitely looks like a bug, thanks for such a detailed report!

[FFroehlich]

Hi Matthias, I am still not sure whether I completely understand the issue you are describing, but I am no longer convinced the specific setup you are describing here is a bug. From my perspective, changing the value of the parameter BW does not constitute a change in volume as this happens outside the scope of a simulation. The same also applies, for example, to the assignment rule for the compartment Vli_plasma, which changes the size of the compartment from 1.0 to Fblood*Vli*(1 - HCT), but does not change the initial concentration of LI__bil_ext to 0.01*Fblood*Vli*(1 - HCT).

However, my understanding is that there are some cases where this is not handled adequately, specifically event mediated changes to concentrations (which we do not support unless there also are rate-rules for the respective compartment, which is probably why we never ran into issues in the sbml testsuite) as well as any petab mediated changes to compartment sizes (which is probably how you ended up looking into these things).

Unfortunately, I don't see an easy way to fix this in AMICI right now and will have to discuss with @dweindl on how this can be best addressed. It's definitely doable, but may require some more substantial changes to the code that are anyways overdue.

[matthiaskoenig]

Hi Fabian.

workaround as a workaround I change the parameters which affect volumes and afterwards set the initial conditions of the affected species again. I.e. I just add the additional changes to the conditions which will reset the values (but for that to work consistently the order of the parameterSpeciesCompartmentIds in the condition table matters:

{
  BW: 80,
  LI__bil_ext: 0.01
}

I.e. I can figure this out programatically for my models and update my conditions accordingly. This solves my main issue, that I have to ensure identical forward simulations for a given condition in the different solvers.

additional info about problem The BW changes the volume via a set of assignment rules.

Vli = BW*FVli * (1-resection_rate)
Vli_plasma = Fblood*Vli*(1 - HCT)

All parameters which via the dependency graph of the assignment rules affect Vli_plasma will have an effect on the concentration of LI__bil_ext via: [LI__bil_ext] = LI__bil_ext/Vli_plasma. I.e. also things such as HCT or Fblood or resection rate will change the LI__bil_ext

The compartments stay constant via the simulations, but after changes to the initial state/parameters the initial compartment volumes must be updated via the assignment rules and initial states depending on the volumes must be updated. I hope this makes things clearer.

By the way, basically what AMICI is doing is what I want, but is just not correct. I.e. I want an actual initial concentration which is 0.01, so I have to enforce this via the correct condition changes:

{
  BW: 80,
  LI__bil_ext: 0.01
}

When I find time I will post a minimal example SBML which shows the issue which should be much easier to debug then this larger model.

[FFroehlich]

Thanks, but the primary issue is that we need different behaviour depending on the format and entries of the condition table, i.e., depending on whether only compartment size or also the concentration is updated. Moreover, the updated initial conditions will depend on the previous size of the volume, which introduces additional parameter (and potentially also state dependencies), which need to be considered to get sensitivities, i.e., gradients right.

All of this needs to be determined during model preprocessing to enable symbolic calculations for sensitivity analysis. The way we implemented condition specific handling of initialization for petab is somewhat of a hack and we need to see whether we extend it or rather go for a proper, a bit more flexible implementation.