assignments_to_experiments: convert values that ramp with time, via a rate rule

petab/v2/converters.py

@@ -674,7 +674,8 @@ def assignments_to_experiments(problem: Problem) -> Problem:
         :data:`DEFAULT_EXPERIMENT_ID`, if needed.
 
     **Scope:** conversion is supported for SBML assignment rules that depend
-    on time only through the conditions of piecewise functions, where
+    on time through the conditions of piecewise functions, and at most
+    linearly through their values, where
 
     * every such condition is an inequality that can be re-arranged
       explicitly into `t </<=/>/>= ...`,
@@ -689,7 +690,13 @@ def assignments_to_experiments(problem: Problem) -> Problem:
     given at different times (`d1 * piecewise(0, t < t1, 1) + d2 *
     piecewise(0, t < t2, 1)`), phases one after another as several pieces or
     nested piecewise functions (`piecewise(b0, t < t1, piecewise(b1, t < t2,
-    b2))`). Each switch that changes the value becomes a period. Anything else
+    b2))`). Each switch that changes the value becomes a period. A value may
+    also change linearly with time within a period (a ramp: `piecewise(r * t
+    + a, t <= t_end, 0)`). A condition change is applied once at the start of
+    its period, so the target of such an assignment becomes a differential
+    entity instead, with the rate rule `d target/dt = <target>_slope`, and
+    each period sets both the target and its slope; an assignment that only
+    switches between constant values stays a constant parameter. Anything else
     raises a :class:`NotImplementedError`.
 
     .. warning::
@@ -774,7 +781,7 @@ def assignments_to_experiments(problem: Problem) -> Problem:
         )
         model_values = {
             s: sbml_model.getParameter(str(s)).getValue()
-            for s in default.free_symbols
+            for s in default.free_symbols - {time}
         }
         # unless a condition sets them, the switching times' symbols keep
         #  their value from the parameter table or else from the model
@@ -786,7 +793,25 @@ def assignments_to_experiments(problem: Problem) -> Problem:
             for s in t_switch.free_symbols
         }
         target = sbml_model.getParameter(target_id)
-        target.setValue(float(default.subs(model_values)))
+        # A value that changes with time within a period (a ramp) cannot be a
+        #  condition change on a constant target, which is applied once at
+        #  the start of its period. The model only makes the target a
+        #  differential entity, `d target/dt = <target>_slope`; the slope's
+        #  value, like the target's, is set by every period that starts one.
+        ramps = expr.xreplace({c: sp.Dummy() for c in switches}).has(time)
+        slope_id = None
+        if ramps:
+            slope_id = f"{target_id}_slope"
+            while sbml_model.getElementBySId(slope_id) is not None:
+                slope_id += "_"
+            slope = sbml_model.createParameter()
+            slope.setId(slope_id)
+            slope.setConstant(False)
+            slope.setValue(0.0)
+            rate_rule = sbml_model.createRateRule()
+            rate_rule.setVariable(target_id)
+            rate_rule.setMath(libsbml.parseL3Formula(slope_id))
+        target.setValue(float(default.subs(time, 0).subs(model_values)))
         target.setConstant(False)
         assignment.removeFromParentAndDelete()
         symbol_ids = {str(s) for s in expr.free_symbols}
@@ -796,7 +821,7 @@ def assignments_to_experiments(problem: Problem) -> Problem:
         for experiment in experiments:
             periods = experiment.sorted_periods
             overrides = {}
-            for i, period in enumerate(periods):
+            for i_period, period in enumerate(periods):
                 changes = problem.get_changes_for_period(period)
                 overrides |= {c.target_id: c.target_value for c in changes}
                 subs = fixed_values | {
@@ -818,14 +843,23 @@ def assignments_to_experiments(problem: Problem) -> Problem:
 
                 if period.is_preequilibration:
                     # the value at time 0, which respects whether a
-                    #  condition holds for a switch at time 0
+                    #  condition holds for a switch at time 0, held there:
+                    #  a ramp does not run during preequilibration
                     value = expr.subs(subs).subs(time, 0)
                     additions.append(
                         (experiment, period.time, target_id, value)
                     )
+                    if ramps:
+                        additions.append(
+                            (experiment, period.time, slope_id, sp.S.Zero)
+                        )
                     continue
 
-                t_end = periods[i + 1].time if i + 1 < len(periods) else sp.oo
+                t_end = (
+                    periods[i_period + 1].time
+                    if i_period + 1 < len(periods)
+                    else sp.oo
+                )
                 start_value = _value_from(
                     expr, switches, time, subs, times, period.time
                 )
@@ -841,9 +875,11 @@ def assignments_to_experiments(problem: Problem) -> Problem:
                 #  the first period, unless that value refers to parameter
                 #  table symbols, for which the model only has a number. The
                 #  switches can become the first periods if they only refer
-                #  to parameter table symbols.
+                #  to parameter table symbols. A ramp's start and slope are
+                #  always set by the tables.
                 redundant = (
-                    i == 0
+                    not ramps
+                    and i_period == 0
                     and (start_value - default).is_zero
                     and not {str(s) for s in default.free_symbols}
                     & parameter_ids
@@ -852,14 +888,38 @@ def assignments_to_experiments(problem: Problem) -> Problem:
                         for _, value in switch_values
                     )
                 )
+                after_preequilibration = (
+                    i_period > 0 and periods[i_period - 1].is_preequilibration
+                )
                 if not redundant and (
-                    i == 0 or {c.target_id for c in changes} & symbol_ids
+                    i_period == 0
+                    or {c.target_id for c in changes} & symbol_ids
+                    # preequilibration held the slope at zero
+                    or (ramps and after_preequilibration)
                 ):
-                    additions.append(
-                        (experiment, period.time, target_id, start_value)
+                    if i_period == 0:
+                        # at the start of an experiment, a change may only
+                        #  refer to parameter table symbols; other symbols
+                        #  are constant model parameters, at their values
+                        start_value = start_value.subs(
+                            {
+                                s: sbml_model.getParameter(str(s)).getValue()
+                                for s in start_value.free_symbols - {time}
+                                if str(s) not in parameter_ids
+                            }
+                        )
+                    additions += _changes(
+                        experiment,
+                        period.time,
+                        target_id,
+                        start_value,
+                        time,
+                        slope_id,
                     )
                 for t, value in switch_values:
-                    additions.append((experiment, t, target_id, value))
+                    additions += _changes(
+                        experiment, t, target_id, value, time, slope_id
+                    )
 
     if converted_ids and any(e.has_preequilibration for e in experiments):
         warnings.warn(
@@ -933,6 +993,16 @@ def assignments_to_experiments(problem: Problem) -> Problem:
     return problem
 
 
+def _changes(experiment, t, target_id, value, time, slope_id):
+    """The changes that start ``value`` at time ``t``: the target set to the
+    value at ``t``, and, for a target that ramps, its slope set to the
+    value's rate of change."""
+    changes = [(experiment, t, target_id, value.subs(time, t))]
+    if slope_id is not None:
+        changes.append((experiment, t, slope_id, sp.diff(value, time)))
+    return changes
+
+
 def _value_from(
     expr: sp.Expr,
     switches: dict,
@@ -993,9 +1063,14 @@ def _supported_time_switches(
                 raise unsupported
             t_switch = sp.solve(condition.gts - condition.lts, time)[0]
             switches[condition] = (t_switch, slope)
-    # time enters only through the conditions
-    if expr.xreplace({c: sp.Dummy() for c in switches}).has(time):
-        raise unsupported
+    # outside the conditions, time enters at most linearly: the value may
+    #  ramp, at a rate that does not itself change with time
+    rest = expr.xreplace({c: sp.Dummy() for c in switches})
+    if rest.has(time) and sp.diff(rest, time).has(time):
+        raise NotImplementedError(
+            "Only values that change at most linearly with time are "
+            f"supported, but the assignment to `{target_id}` is `{expr}`."
+        )
     # no period can start at an estimated time
     if estimated := sorted(
         {

tests/v2/test_conversion.py

@@ -1,6 +1,7 @@
 import logging
 from pathlib import Path
 
+import libsbml
 import numpy as np
 import pandas as pd
 import pytest
@@ -939,7 +940,8 @@ def test_assignments_to_experiments_preequilibration():
         ("piecewise(0, time == 10, 1)", "single"),
         ("piecewise(0, time < 10 && time > 5, 1)", "single"),
         ("piecewise(0, time^2 < 10, 1)", "single"),
-        ("piecewise(time, time < 10, 1)", "single"),
+        ("piecewise(time^2, time < 10, 1)", "linearly"),
+        ("piecewise(exp(time), time < 10, 1)", "linearly"),
         (
             "piecewise(0, time < t_dose, 1)",
             r"estimated parameters \['t_dose'\]",
@@ -1178,6 +1180,149 @@ def test_assignments_to_experiments_condition_specific_doses():
     }
 
 
+def _rate_rule(problem, target_id):
+    rule = problem.model.sbml_model.getRateRuleByVariable(target_id)
+    return None if rule is None else libsbml.formulaToL3String(rule.getMath())
+
+
+def _changes_by_period(problem):
+    """{(experiment ID, time): {target ID: target value}}"""
+    conditions = problem.condition_df.set_index(v2.C.CONDITION_ID)
+    out = {}
+    for _, row in problem.experiment_df.iterrows():
+        change = conditions.loc[[row[v2.C.CONDITION_ID]]]
+        out.setdefault((row[v2.C.EXPERIMENT_ID], row[v2.C.TIME]), {}).update(
+            {
+                r[v2.C.TARGET_ID]: float(r[v2.C.TARGET_VALUE])
+                for _, r in change.iterrows()
+            }
+        )
+    return out
+
+
+def test_assignments_to_experiments_ramp():
+    """Test that a value changing linearly with time becomes a rate rule
+    whose slope each period sets."""
+    problem = Problem()
+    ant_model = """
+        model conversion
+        species A, B;
+        A = 10;
+        B = 0;
+        r = 0.5; a = 1;
+        k1 := piecewise(r * time + a, time <= 10, 0);
+        k2 = 0.5;
+        R1: A -> B; k1 * A;
+        R2: B -> A; k2 * B;
+        end
+    """
+    problem.model = SbmlModel.from_antimony(ant_model)
+    problem.add_observable("obs_a", formula="A", noise_formula=1.0)
+    problem.add_measurement("obs_a", time=20.0, measurement=5.0)
+    problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3)
+
+    converted = assignments_to_experiments(problem)
+
+    sbml = converted.model.sbml_model
+    # the model only makes the target differential; the slope is a
+    #  placeholder until the tables set it
+    assert _rate_rule(converted, "k1") == "k1_slope"
+    assert sbml.getParameter("k1_slope").getValue() == 0.0
+    # the ramp starts at time 0, and at the switch, the target and its slope
+    #  drop to zero
+    assert _changes_by_period(converted) == {
+        (DEFAULT_EXPERIMENT_ID, 0.0): {"k1": 1.0, "k1_slope": 0.5},
+        (DEFAULT_EXPERIMENT_ID, 10.0): {"k1": 0.0, "k1_slope": 0.0},
+    }
+
+
+def test_assignments_to_experiments_condition_specific_ramps():
+    """Test pulses and ramps whose start, rate and end are set by the
+    condition table, as in Fujita_SciSignal2010."""
+    problem = Problem()
+    ant_model = """
+        model conversion
+        species A, B;
+        A = 10;
+        B = 0;
+        EGF_0 = 0; EGF_rate = 0; EGF_end = 0;
+        k1 := piecewise(EGF_rate * time + EGF_0, time <= EGF_end, 0);
+        k2 = 0.5;
+        R1: A -> B; k1 * A;
+        R2: B -> A; k2 * B;
+        end
+    """
+    problem.model = SbmlModel.from_antimony(ant_model)
+    # a pulse: a constant dose withdrawn at t = 60
+    problem.add_condition("pulse", EGF_0=1.0, EGF_rate=0.0, EGF_end=60.0)
+    # a ramp: growing from zero at 0.01 per unit time until t = 3600
+    problem.add_condition("ramp", EGF_0=0.0, EGF_rate=0.01, EGF_end=3600.0)
+    for c in ("pulse", "ramp"):
+        problem.add_experiment(f"experiment____{c}", 0, c)
+        problem.add_measurement(
+            "obs_a",
+            experiment_id=f"experiment____{c}",
+            time=30.0,
+            measurement=5.0,
+        )
+    problem.add_observable("obs_a", formula="A", noise_formula=1.0)
+    problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3)
+
+    converted = assignments_to_experiments(problem)
+
+    assert _rate_rule(converted, "k1") == "k1_slope"
+    assert _changes_by_period(converted) == {
+        ("experiment____pulse", 0.0): {"k1": 1.0, "k1_slope": 0.0},
+        ("experiment____pulse", 60.0): {"k1": 0.0, "k1_slope": 0.0},
+        ("experiment____ramp", 0.0): {"k1": 0.0, "k1_slope": 0.01},
+        ("experiment____ramp", 3600.0): {"k1": 0.0, "k1_slope": 0.0},
+    }
+
+
+def test_assignments_to_experiments_ramp_preequilibration():
+    """Test that a ramp is held at its value at time 0 during
+    preequilibration and starts at time 0."""
+    problem = Problem()
+    ant_model = """
+        model conversion
+        species A, B;
+        A = 10;
+        B = 0;
+        r = 0;
+        k1 := piecewise(r * time + 1, time <= 10, 0);
+        k2 = 0.5;
+        R1: A -> B; k1 * A;
+        R2: B -> A; k2 * B;
+        end
+    """
+    problem.model = SbmlModel.from_antimony(ant_model)
+    problem.add_condition("pre", r=0.0)
+    problem.add_condition("main", r=2.0)
+    problem.add_experiment(
+        "experiment__pre___main", v2.C.TIME_PREEQUILIBRATION, "pre", 0, "main"
+    )
+    problem.add_observable("obs_a", formula="A", noise_formula=1.0)
+    problem.add_measurement(
+        "obs_a",
+        experiment_id="experiment__pre___main",
+        time=20.0,
+        measurement=5.0,
+    )
+    problem.add_parameter("k2", estimate=True, lb=1e-3, ub=1e3)
+
+    with pytest.warns(UserWarning, match="during preequilibration"):
+        converted = assignments_to_experiments(problem)
+
+    assert _changes_by_period(converted) == {
+        ("experiment__pre___main", v2.C.TIME_PREEQUILIBRATION): {
+            "k1": 1.0,
+            "k1_slope": 0.0,
+        },
+        ("experiment__pre___main", 0.0): {"k1": 1.0, "k1_slope": 2.0},
+        ("experiment__pre___main", 10.0): {"k1": 0.0, "k1_slope": 0.0},
+    }
+
+
 @pytest.mark.parametrize("t_switch", [0, 10])
 @pytest.mark.parametrize(
     "condition, holds_at_switch, holds_before",