Generation of condition-specific SBML models (#108)

New function `petab.get_model_for_condition` to generate condition-specific SBML models based on a PEtab problem.

Co-authored-by: Dilan Pathirana <[email protected]>

Author: dweindl

petab/sbml.py

@@ -1,20 +1,24 @@
 """Functions for interacting with SBML models"""
 import logging
 import re
+from numbers import Number
 from pathlib import Path
-from typing import Any, Dict, List, Tuple, Union
+from typing import Any, Dict, List, Optional, Tuple, Union
 from warnings import warn
 
 import libsbml
 from pandas.io.common import get_handle, is_file_like, is_url
 
+import petab
+
 logger = logging.getLogger(__name__)
 __all__ = ['add_global_parameter',
            'add_model_output',
            'add_model_output_sigma',
            'add_model_output_with_sigma',
            'assignment_rules_to_dict',
            'create_assigment_rule',
+           'get_model_for_condition',
            'get_model_parameters',
            'get_observables',
            'get_sbml_model',
@@ -481,3 +485,118 @@ def load_sbml_from_file(
     sbml_model = sbml_document.getModel()
 
     return sbml_reader, sbml_document, sbml_model
+
+
+def get_model_for_condition(
+        petab_problem: "petab.Problem",
+        sim_condition_id: str = None,
+        preeq_condition_id: Optional[str] = None,
+) -> Tuple[libsbml.SBMLDocument, libsbml.Model]:
+    """Create an SBML model for the given condition.
+
+    Creates a copy of the model and updates parameters according to the PEtab
+    files. Estimated parameters are set to their ``nominalValue``.
+    Observables defined in the observables table are not added to the model.
+
+    :param petab_problem: PEtab problem
+    :param sim_condition_id: Simulation ``conditionId`` for which to generate a
+        model
+    :param preeq_condition_id: Preequilibration ``conditionId`` of the settings
+        for which to generate a model. This is only used to determine the
+        relevant output parameter overrides. Preequilibration is not encoded
+        in the resulting model.
+    :return: The generated SBML document, and SBML model
+    """
+    condition_dict = {petab.SIMULATION_CONDITION_ID: sim_condition_id}
+    if preeq_condition_id:
+        condition_dict[petab.PREEQUILIBRATION_CONDITION_ID] = \
+            preeq_condition_id
+    cur_measurement_df = petab.measurements.get_rows_for_condition(
+        measurement_df=petab_problem.measurement_df,
+        condition=condition_dict,
+    )
+    parameter_map, scale_map = \
+        petab.parameter_mapping.get_parameter_mapping_for_condition(
+            condition_id=sim_condition_id,
+            is_preeq=False,
+            cur_measurement_df=cur_measurement_df,
+            sbml_model=petab_problem.sbml_model,
+            condition_df=petab_problem.condition_df,
+            parameter_df=petab_problem.parameter_df,
+            warn_unmapped=True,
+            scaled_parameters=False,
+            fill_fixed_parameters=True,
+            # will only become problematic once the observable and noise terms
+            #  are added to the model
+            allow_timepoint_specific_numeric_noise_parameters=True,
+        )
+    # create a copy of the model
+    sbml_doc = petab_problem.sbml_model.getSBMLDocument().clone()
+    sbml_model = sbml_doc.getModel()
+
+    # fill in parameters
+    def get_param_value(parameter_id: str):
+        """Parameter value from mapping or nominal value"""
+        mapped_value = parameter_map.get(parameter_id)
+        if not isinstance(mapped_value, str):
+            return mapped_value
+
+        # estimated parameter, look up in nominal parameters
+        return petab_problem.parameter_df.loc[mapped_value,
+                                              petab.NOMINAL_VALUE]
+
+    def remove_rules(target_id: str):
+        if sbml_model.removeRuleByVariable(target_id):
+            warn("An SBML rule was removed to set the component "
+                 f"{target_id} to a constant value.")
+        sbml_model.removeInitialAssignment(target_id)
+
+    for parameter in sbml_model.getListOfParameters():
+        new_value = get_param_value(parameter.getId())
+        if new_value:
+            parameter.setValue(new_value)
+            # remove rules that would override that value
+            remove_rules(parameter.getId())
+
+    # set concentrations for any overridden species
+    for component_id in petab_problem.condition_df:
+        sbml_species = sbml_model.getSpecies(component_id)
+        if not sbml_species:
+            continue
+
+        # remove any rules overriding that species' initials
+        remove_rules(component_id)
+
+        # set initial concentration/amount
+        new_value = petab.to_float_if_float(
+            petab_problem.condition_df.loc[sim_condition_id, component_id])
+        if not isinstance(new_value, Number):
+            # parameter reference in condition table
+            new_value = get_param_value(new_value)
+
+        if sbml_species.isSetInitialAmount() \
+            or (sbml_species.getHasOnlySubstanceUnits()
+                and not sbml_species.isSetInitialConcentration()):
+            sbml_species.setInitialAmount(new_value)
+        else:
+            sbml_species.setInitialConcentration(new_value)
+
+    # set compartment size for any compartments in the condition table
+    for component_id in petab_problem.condition_df:
+        sbml_compartment = sbml_model.getCompartment(component_id)
+        if not sbml_compartment:
+            continue
+
+        # remove any rules overriding that compartment's size
+        remove_rules(component_id)
+
+        # set initial concentration/amount
+        new_value = petab.to_float_if_float(
+            petab_problem.condition_df.loc[sim_condition_id, component_id])
+        if not isinstance(new_value, Number):
+            # parameter reference in condition table
+            new_value = get_param_value(new_value)
+
+        sbml_compartment.setSize(new_value)
+
+    return sbml_doc, sbml_model

tests/test_sbml.py

@@ -2,6 +2,8 @@
 import sys
 import os
 
+import pandas as pd
+
 sys.path.append(os.getcwd())
 import petab  # noqa: E402
 
@@ -31,3 +33,83 @@ def test_assignment_rules_to_dict():
     assert actual == expected
     assert model.getAssignmentRuleByVariable('observable_1') is None
     assert len(model.getListOfParameters()) == 0
+
+
+def test_get_condition_specific_models():
+    """Test for petab.sbml.get_condition_specific_models"""
+    # Create test model and data files
+    sbml_document = libsbml.SBMLDocument(3, 1)
+    sbml_model = sbml_document.createModel()
+
+    c = sbml_model.createCompartment()
+    c.setId("compartment_1")
+    c.setSize(1)
+
+    for i in range(1, 4):
+        p = sbml_model.createParameter()
+        p.setId(f"parameter_{i}")
+        p.setValue(i)
+
+    for i in range(1, 4):
+        s = sbml_model.createSpecies()
+        s.setId(f"species_{i}")
+        s.setInitialConcentration(10 * i)
+
+    ia = sbml_model.createInitialAssignment()
+    ia.setSymbol("species_2")
+    ia.setMath(libsbml.parseL3Formula("25"))
+
+    condition_df = pd.DataFrame({
+        petab.CONDITION_ID: ["condition_1"],
+        "parameter_3": ['parameter_2'],
+        "species_1": [15],
+        "species_2": [25],
+        "species_3": ['parameter_1'],
+        "compartment_1": [2],
+    })
+    condition_df.set_index([petab.CONDITION_ID], inplace=True)
+
+    observable_df = pd.DataFrame({
+        petab.OBSERVABLE_ID: ["observable_1"],
+        petab.OBSERVABLE_FORMULA: ["2 * species_1"],
+    })
+    observable_df.set_index([petab.OBSERVABLE_ID], inplace=True)
+
+    measurement_df = pd.DataFrame({
+        petab.OBSERVABLE_ID: ["observable_1"],
+        petab.SIMULATION_CONDITION_ID: ["condition_1"],
+        petab.TIME: [0.0],
+    })
+
+    parameter_df = pd.DataFrame({
+        petab.PARAMETER_ID: ["parameter_1", "parameter_2"],
+        petab.PARAMETER_SCALE: [petab.LOG10] * 2,
+        petab.NOMINAL_VALUE: [1.25, 2.25],
+        petab.ESTIMATE: [0, 1],
+    })
+    parameter_df.set_index([petab.PARAMETER_ID], inplace=True)
+
+    petab_problem = petab.Problem(
+        sbml_model=sbml_model,
+        condition_df=condition_df,
+        observable_df=observable_df,
+        measurement_df=measurement_df,
+        parameter_df=parameter_df
+    )
+
+    # Actual test
+    condition_doc, condition_model = petab.get_model_for_condition(
+        petab_problem, "condition_1")
+
+    assert condition_model.getSpecies(
+        "species_1").getInitialConcentration() == 15
+    assert condition_model.getSpecies(
+        "species_2").getInitialConcentration() == 25
+    assert condition_model.getSpecies(
+        "species_3").getInitialConcentration() == 1.25
+    assert len(condition_model.getListOfInitialAssignments()) == 0, \
+        "InitialAssignment not removed"
+    assert condition_model.getCompartment("compartment_1").getSize() == 2.0
+    assert condition_model.getParameter("parameter_1").getValue() == 1.25
+    assert condition_model.getParameter("parameter_2").getValue() == 2.25
+    assert condition_model.getParameter("parameter_3").getValue() == 2.25