Update to changed change semantics

Author: dweindl

petab/v2/converters.py

@@ -6,6 +6,7 @@
 from copy import deepcopy
 
 import libsbml
+import sympy as sp
 from sbmlmath import sbml_math_to_sympy, set_math
 
 from .core import (
@@ -31,6 +32,8 @@ class ExperimentsToEventsConverter:
     If the model already contains events, PEtab events are added with a higher
     priority than the existing events to guarantee that PEtab condition changes
     are applied before any pre-existing assignments.
+    This requires that all event priorities in the original model are numeric
+    constants.
 
     The PEtab problem must not contain any identifiers starting with
     ``_petab``.
@@ -92,9 +95,8 @@ def __init__(self, problem: Problem, default_priority: float = None):
         self._default_priority = default_priority
         self._preprocess()
 
-    def _get_experiment_indicator_condition_id(
-        self, experiment_id: str
-    ) -> str:
+    @staticmethod
+    def _get_experiment_indicator_condition_id(experiment_id: str) -> str:
         """Get the condition ID for the experiment indicator parameter."""
         return f"_petab_experiment_condition_{experiment_id}"
 
@@ -214,7 +216,7 @@ def _convert_experiment(self, experiment: Experiment) -> None:
             )
         add_sbml_parameter(model, id_=exp_ind_id, constant=False, value=0)
         kept_periods = []
-        for i_period, period in enumerate(experiment.periods):
+        for i_period, period in enumerate(experiment.sorted_periods):
             if period.is_preequilibration:
                 # pre-equilibration cannot be represented in SBML,
                 #  so we need to keep this period in the Problem.
@@ -229,17 +231,20 @@ def _convert_experiment(self, experiment: Experiment) -> None:
                 #  or the only non-equilibration period (handled above)
                 continue
 
+            # TODO: If this is the first period, we may want to implement
+            #  the changes as initialAssignments instead of events.
+            #  That way, we don't need to worry about any event priorities
+            #  and tools that don't support event priorities can still handle
+            #  single-period experiments.
+
             ev = self._create_period_start_event(
                 experiment=experiment,
                 i_period=i_period,
                 period=period,
             )
             self._create_event_assignments_for_period(
                 ev,
-                [
-                    self._new_problem[condition_id]
-                    for condition_id in period.condition_ids
-                ],
+                self._new_problem.get_changes_for_period(period),
             )
 
         if len(kept_periods) > 2:
@@ -326,33 +331,120 @@ def get_experiment_indicator(experiment_id: str) -> str:
 
     @staticmethod
     def _create_event_assignments_for_period(
-        event: libsbml.Event, conditions: list[Condition]
+        event: libsbml.Event, changes: list[Change]
     ) -> None:
-        """Create an event assignments for a given period."""
-        for condition in conditions:
-            for change in condition.changes:
-                ExperimentsToEventsConverter._change_to_event_assignment(
-                    change, event
+        """Create event assignments for a given period.
+
+        Converts PEtab ``Change``s to equivalent SBML event assignments.
+
+        Note that the SBML event assignment formula is not necessarily the same
+        as the `targetValue` in PEtab.
+        In SBML, concentrations are treated as derived quantities.
+        Therefore, changing the size of a compartment will update the
+        concentrations of all contained concentration-based species.
+        In PEtab, such a change would not automatically update the species
+        concentrations, but only the compartment size.
+
+        Therefore, to correctly implement a PEtab change of a compartment size
+        in SBML, we need to compensate for the automatic update of species
+        concentrations by adding event assignments for all contained
+        concentration-based species.
+
+        :param event: The SBML event to which the assignments should be added.
+        :param changes: The PEtab condition changes that are to be applied
+            at the start of the period.
+        """
+        _add_assignment = ExperimentsToEventsConverter._add_assignment
+        sbml_model = event.getModel()
+        # collect IDs of compartments that are changed in this period
+        changed_compartments = {
+            change.target_id
+            for change in changes
+            if sbml_model.getElementBySId(change.target_id) is not None
+            and sbml_model.getElementBySId(change.target_id).getTypeCode()
+            == libsbml.SBML_COMPARTMENT
+        }
+
+        for change in changes:
+            sbml_target = sbml_model.getElementBySId(change.target_id)
+
+            if sbml_target is None:
+                raise ValueError(
+                    f"Cannot create event assignment for change of "
+                    f"`{change.target_id}`: No such entity in the SBML model."
                 )
 
+            target_type = sbml_target.getTypeCode()
+            if target_type == libsbml.SBML_COMPARTMENT:
+                # handle the actual compartment size change
+                _add_assignment(event, change.target_id, change.target_value)
+
+                # Changing a compartment size affects all contained
+                #  concentration-based species - we need to add event
+                #  assignments for those to compensate for the automatic
+                #  update of their concentrations.
+                # The event assignment will set the concentration to
+                #   new_conc = assigned_amount / new_volume
+                #            = assigned_conc * old_volume / new_volume
+                #   <=> assigned_conc = new_conc * new_volume / old_volume
+                # Therefore, the event assignment is not just `new_conc`,
+                #  but `new_conc * new_volume / old_volume`.
+
+                # concentration-based species in the changed compartment
+                conc_species = [
+                    species.getId()
+                    for species in sbml_model.getListOfSpecies()
+                    if species.getCompartment() == change.target_id
+                    and not species.getHasOnlySubstanceUnits()
+                ]
+                for species_id in conc_species:
+                    if species_change := next(
+                        (c for c in changes if c.target_id == species_id), None
+                    ):
+                        # there is an explicit change for this species
+                        #  in this period
+                        new_conc = species_change.target_value
+                    else:
+                        # no explicit change, use the pre-event concentration
+                        new_conc = sp.Symbol(species_id)
+
+                    _add_assignment(
+                        event,
+                        species_id,
+                        # new_conc * new_volume / old_volume
+                        new_conc
+                        * change.target_value
+                        / sp.Symbol(change.target_id),
+                    )
+            elif (
+                target_type != libsbml.SBML_SPECIES
+                or sbml_target.getCompartment() not in changed_compartments
+                or sbml_target.getHasOnlySubstanceUnits() is True
+            ):
+                # Handle any changes other than compartments and
+                #  concentration-based species inside resized compartments
+                #  that we were already handled above.
+                # Those translate directly to event assignments.
+                _add_assignment(event, change.target_id, change.target_value)
+
     @staticmethod
-    def _change_to_event_assignment(
-        change: Change, event: libsbml.Event
+    def _add_assignment(
+        event: libsbml.Event, target_id: str, target_value: sp.Basic
     ) -> None:
-        """Convert a PEtab ``Change``  to an SBML event assignment."""
+        """Add a single event assignment to the given event
+        and apply any necessary changes to the model."""
         sbml_model = event.getModel()
-
         ea = event.createEventAssignment()
-        ea.setVariable(change.target_id)
-        set_math(ea, change.target_value)
+        ea.setVariable(target_id)
+        set_math(ea, target_value)
 
         # target needs const=False, and target may not exist yet
         #  (e.g., in case of output parameters added in the observable
         #  table)
-        target = sbml_model.getElementBySId(change.target_id)
+        target = sbml_model.getElementBySId(target_id)
         if target is None:
             add_sbml_parameter(
-                sbml_model, id_=change.target_id, constant=False, value=0
+                sbml_model, id_=target_id, constant=False, value=0
             )
         else:
             # We can safely change the `constant` attribute of the target.
@@ -362,7 +454,7 @@ def _change_to_event_assignment(
         # the target value may depend on parameters that are only
         #  introduced in the PEtab parameter table - those need
         #  to be added to the model
-        for sym in change.target_value.free_symbols:
+        for sym in target_value.free_symbols:
             if sbml_model.getElementBySId(sym.name) is None:
                 add_sbml_parameter(
                     sbml_model, id_=sym.name, constant=True, value=0

pyproject.toml

@@ -37,6 +37,7 @@ maintainers = [
 [project.optional-dependencies]
 tests = [
     "antimony>=2.14.0",
+    "copasi-basico>=0.85",
     "pysb",
     "pytest",
     "pytest-cov",

tests/v2/test_converters.py

@@ -1,5 +1,7 @@
 from math import inf
 
+import pandas as pd
+
 from petab.v2 import Change, Condition, Experiment, ExperimentPeriod, Problem
 from petab.v2.converters import ExperimentsToEventsConverter
 from petab.v2.models.sbml_model import SbmlModel
@@ -74,3 +76,156 @@ def test_experiments_to_events_converter():
             ],
         ),
     ]
+
+
+def test_simulate_experiment_to_events():
+    """
+    Convert PEtab experiment to SBML events and compare BasiCO simulation
+    results.
+    """
+    import basico
+
+    # the basic model for the PEtab problem
+    ant_model1 = """
+    compartment comp1 = 10
+    compartment comp2 = 2
+    # concentration-based species
+    species s1c_comp1 in comp1 = 1
+    species s1c_comp2 in comp2 = 2
+    species s2c_comp1 in comp1 = 3
+    species s2c_comp2 in comp2 = 4
+    # amount-based species
+    # (note: the initial values are concentrations nonetheless)
+    substanceOnly species s3a_comp1 in comp1 = 5
+    substanceOnly species s3a_comp2 in comp2 = 6
+    substanceOnly species s4a_comp1 in comp1 = 7
+    substanceOnly species s4a_comp2 in comp2 = 8
+
+    # something dynamic
+    some_species in comp1 = 0
+    some_species' = 1
+
+    # set time-derivatives, otherwise BasiCO won't include them in the result
+    s1c_comp1' = 0
+    s1c_comp2' = 0
+    s2c_comp1' = 0
+    s2c_comp2' = 0
+    s3a_comp1' = 0
+    s3a_comp2' = 0
+    s4a_comp1' = 0
+    s4a_comp2' = 0
+    """
+
+    # append events, equivalent to the expected PEtab conversion result
+    ant_model_expected = (
+        ant_model1
+        + """
+    # resize compartment
+    # The size of comp1 should be set to 20, the concentrations of the
+    # contained concentration-based species and the amounts of the amount-based
+    # species should remain unchanged. comp2 and everything therein is
+    # unaffected.
+    # I.e., post-event:
+    #   s1c_comp1 = 1, s2c_comp1 = 3, s3a_comp1 = 5, s4a_comp1 = 7
+    at time >= 1:
+        comp1 = 20,
+        s1c_comp1 = s1c_comp1 * 20 / comp1,
+        s2c_comp1 = s2c_comp1 * 20 / comp1;
+
+    # resize compartment *and* reassign concentration
+    # The size of comp2 should be set to 4, the concentration/amount of
+    # s1c_comp2/s3a_comp2 should be set to the given values,
+    # the amounts for amount-based and concentrations for concentration-based
+    # other species in comp2 should remain unchanged.
+    # I.e., post-event:
+    #   comp2 = 4
+    #   s1c_comp2 = 5, s3a_comp2 = 16,
+    #   s2c_comp2 = 4 (unchanged), s4a_comp2 = 8 (unchanged)
+    # The post-event concentrations of concentration-based species are
+    # (per SBML):
+    #   new_conc = assigned_amount / new_volume
+    #            = assigned_conc * old_volume / new_volume
+    #   <=> assigned_conc = new_conc * new_volume / old_volume
+    # The post-event amounts of amount-based species are:
+    #   new_amount = assigned_amount (independent of volume change)
+    at time >= 5:
+        comp2 = 4,
+        s3a_comp2 = 16,
+        s1c_comp2 = 5 * 4 / comp2,
+        s2c_comp2 = s2c_comp2 * 4 / comp2;
+    """
+    )
+
+    # simulate expected model in BasiCO
+    sbml_expected = SbmlModel.from_antimony(ant_model_expected).to_sbml_str()
+    basico.load_model(sbml_expected)
+    # output timepoints (initial, pre-/post-event, ...)
+    timepoints = [0, 0.9, 1.1, 4.9, 5.1, 10]
+    # Simulation will return all species as concentrations
+    df_expected = basico.run_time_course(values=timepoints)
+    # fmt: off
+    assert (
+            df_expected
+            == pd.DataFrame(
+            {'Values[some_species]': {0.0: 0.0, 0.9: 0.9,
+                                      1.1: 1.0999999999999996, 4.9: 4.9,
+                                      5.1: 5.100000000000001, 10.0: 10.0},
+             's1c_comp1': {0.0: 1.0, 0.9: 1.0, 1.1: 1.0, 4.9: 1.0, 5.1: 1.0,
+                           10.0: 1.0},
+             's2c_comp1': {0.0: 3.0, 0.9: 3.0, 1.1: 3.0, 4.9: 3.0, 5.1: 3.0,
+                           10.0: 3.0},
+             's3a_comp1': {0.0: 5.0, 0.9: 5.0, 1.1: 2.5, 4.9: 2.5, 5.1: 2.5,
+                           10.0: 2.5},
+             's4a_comp1': {0.0: 7.0, 0.9: 7.0, 1.1: 3.5, 4.9: 3.5, 5.1: 3.5,
+                           10.0: 3.5},
+             's1c_comp2': {0.0: 2.0, 0.9: 2.0, 1.1: 2.0, 4.9: 2.0, 5.1: 5.0,
+                           10.0: 5.0},
+             's2c_comp2': {0.0: 4.0, 0.9: 4.0, 1.1: 4.0, 4.9: 4.0, 5.1: 4.0,
+                           10.0: 4.0},
+             's3a_comp2': {0.0: 6.0, 0.9: 6.0, 1.1: 6.0, 4.9: 6.0, 5.1: 4.0,
+                           10.0: 4.0},
+             's4a_comp2': {0.0: 8.0, 0.9: 8.0, 1.1: 8.0, 4.9: 8.0, 5.1: 4.0,
+                           10.0: 4.0},
+             'Compartments[comp1]': {0.0: 10.0, 0.9: 10.0, 1.1: 20.0,
+                                     4.9: 20.0, 5.1: 20.0, 10.0: 20.0},
+             'Compartments[comp2]': {0.0: 2.0, 0.9: 2.0, 1.1: 2.0, 4.9: 2.0,
+                                     5.1: 4.0, 10.0: 4.0}}
+        )
+    ).all().all()
+    # fmt: on
+
+    # construct PEtab test problem
+    problem = Problem()
+    problem.model = SbmlModel.from_antimony(ant_model1)
+    problem.add_condition("c1", comp1=20)
+    problem.add_condition("c2", comp2=4, s1c_comp2=5, s3a_comp2=16)
+    problem.add_experiment("e1", 1, "c1", 5, "c2")
+    problem.assert_valid()
+
+    # convert PEtab experiments to SBML events and simulate in BasiCO
+    converter = ExperimentsToEventsConverter(problem)
+    converted = converter.convert()
+    # set experiment indicator to simulate experiment "e1"
+    converted.model.sbml_model.getParameter(
+        "_petab_experiment_indicator_e1"
+    ).setValue(1)
+    sbml_actual = converted.model.to_sbml_str()
+    basico.load_model(sbml_actual)
+    df_actual = basico.run_time_course(values=timepoints)
+
+    # compare results
+    with pd.option_context(
+        "display.max_rows",
+        None,
+        "display.max_columns",
+        None,
+        "display.width",
+        None,
+    ):
+        print("Expected:")
+        print(df_expected)
+        print("Actual:")
+        print(df_actual)
+
+    for col in df_expected.columns:
+        assert (df_expected[col] == df_actual[col]).all()