diff --git a/petab/v2/converters.py b/petab/v2/converters.py index 8dc710f2..6b2d8217 100644 --- a/petab/v2/converters.py +++ b/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 = _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 = _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( { diff --git a/tests/v2/test_conversion.py b/tests/v2/test_conversion.py index 8cdc5237..6f621bb5 100644 --- a/tests/v2/test_conversion.py +++ b/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",