Add Blasi model selection problem for motif-specific parameters (#157)

Author: dilpath

Benchmark-Models/Blasi_CellSystems2016/petab_select/01_create_output_dirs.sh

@@ -0,0 +1,3 @@
+mkdir -p output/model
+mkdir output/petab
+mkdir output/select

Benchmark-Models/Blasi_CellSystems2016/petab_select/02_create_sbml.py

@@ -0,0 +1,66 @@
+import libsbml
+
+from model_info import (
+    all_species_dict,
+    reactions_dict,
+    parameters_dict,
+)
+
+
+# Create model
+sbml_document = libsbml.SBMLDocument(3, 2)
+sbml_model = sbml_document.createModel()
+sbml_model.setId('Blasi_CellSystems2016__select_motif_specific')
+sbml_model.setName('Blasi_CellSystems2016__select_motif_specific')
+
+# Add compartment
+compartment_id = 'compartment'
+compartment = sbml_model.createCompartment()
+compartment.setId(compartment_id)
+compartment.setName(compartment_id)
+compartment.setSize(1)
+compartment.setConstant(True)
+
+# Add species
+for species_dict in all_species_dict.values():
+    species = sbml_model.createSpecies()
+    species.setId(species_dict['id'])
+    species.setName(species_dict['name'])
+    species.setInitialConcentration(species_dict['initialConcentration'])
+    species.setCompartment(compartment_id)
+    species.setConstant(False)
+    species.setBoundaryCondition(False)
+    species.setHasOnlySubstanceUnits(False)
+
+# Add parameters
+for parameter_dict in parameters_dict.values():
+    parameter = sbml_model.createParameter()
+    parameter.setId(parameter_dict['id'])
+    parameter.setName(parameter_dict['name'])
+    parameter.setValue(parameter_dict['value'])
+    parameter.setConstant(True)
+
+# Add reactions
+for reaction_dict in reactions_dict.values():
+    # Create reaction
+    reaction = sbml_model.createReaction()
+    reaction.setId(reaction_dict['id'])
+    reaction.setName(reaction_dict['name'])
+    reaction.setReversible(False)
+
+    # Add reactant
+    reactant = reaction.createReactant()
+    reactant.setSpecies(reaction_dict['reactant_id'])
+    reactant.setConstant(True)
+
+    # Add product
+    product = reaction.createProduct()
+    product.setSpecies(reaction_dict['product_id'])
+    product.setConstant(True)
+
+    # Add kinetic
+    math_ast = libsbml.parseL3Formula(reaction_dict['formula'])
+    kinetic_law = reaction.createKineticLaw()
+    kinetic_law.setMath(math_ast)
+
+libsbml.writeSBMLToFile(sbml_document, 'output/model/model.xml')

Benchmark-Models/Blasi_CellSystems2016/petab_select/03_create_petab.py

@@ -0,0 +1,159 @@
+import pandas as pd
+import petab
+import shutil
+
+from petab import (
+    CONDITION_ID,
+    OBSERVABLE_ID,
+    MEASUREMENT,
+    OBSERVABLE_NAME,
+    OBSERVABLE_FORMULA,
+    NOISE_FORMULA,
+    OBSERVABLE_TRANSFORMATION,
+    NOISE_DISTRIBUTION,
+    LOG,
+    LOG10,
+    NORMAL,
+    SIMULATION_CONDITION_ID,
+    TIME,
+    PARAMETER_ID,
+    PARAMETER_NAME,
+    LIN,
+    PARAMETER_SCALE,
+    LOWER_BOUND,
+    UPPER_BOUND,
+    NOMINAL_VALUE,
+    ESTIMATE,
+)
+
+
+from model_info import (
+    get_ac_state,
+    get_motif,
+    get_species_id,
+    lysines,
+    motif_specific_parameters,
+    # switch_parameters,
+    common_parameters,
+)
+
+
+petab_problem0 = petab.Problem.from_yaml('../Blasi_CellSystems2016.yaml')
+measurements0 = dict(zip(petab_problem0.measurement_df[OBSERVABLE_ID], petab_problem0.measurement_df[MEASUREMENT]))
+
+condition_id = 'condition'
+noise = 'sigma'
+
+
+def convert_observable_id(id0: str) -> str:
+    motif0 = id0[len('observable_'):]
+    if motif0 == '0ac':
+        ac_state = tuple()
+    elif motif0 == '4ac':
+        ac_state = get_ac_state(lysines)
+    else:
+        ac_state = get_ac_state([f'k{i:0>2}' for i in motif0.split('k') if i])
+    species_id = get_species_id(ac_state)
+    observable_id = 'observable_' + species_id
+    observable_dict = {
+        OBSERVABLE_ID: 'observable_' + species_id,
+        OBSERVABLE_NAME: 'y_{' + get_motif(ac_state) + '}',
+        OBSERVABLE_FORMULA: species_id,
+        OBSERVABLE_TRANSFORMATION: LOG,
+        NOISE_FORMULA: noise,
+        NOISE_DISTRIBUTION: NORMAL,
+    }
+    return observable_dict
+
+
+# Create observables
+observable_dicts = []
+observable_id_mapping = {}
+observable_ids0 = sorted(set(measurements0))
+for observable_id0 in observable_ids0:
+    observable_dict = convert_observable_id(observable_id0)
+    observable_id_mapping[observable_id0] = observable_dict[OBSERVABLE_ID]
+    observable_dicts.append(observable_dict)
+
+# Create measurements
+measurement_dicts = []
+for observable_id0, value0 in measurements0.items():
+    observable_id = observable_id_mapping[observable_id0]
+    measurement_dict = {
+        OBSERVABLE_ID: observable_id,
+        SIMULATION_CONDITION_ID: condition_id,
+        TIME: 'inf',
+        MEASUREMENT: value0,
+    }
+    measurement_dicts.append(measurement_dict)
+
+# Create parameters
+parameter_dicts = []
+## Motif-specific
+for parameter_dict0 in motif_specific_parameters.values():
+    parameter_dict = {
+        PARAMETER_ID: parameter_dict0['id'],
+        PARAMETER_NAME: parameter_dict0['name'],
+        PARAMETER_SCALE: LOG10,
+        LOWER_BOUND: '1e-12',
+        UPPER_BOUND: '1e3',
+        NOMINAL_VALUE: 1,
+        ESTIMATE: 0,
+    }
+    parameter_dicts.append(parameter_dict)
+# ## Switch
+# for parameter_dict0 in switch_parameters.values():
+#     parameter_dict = {
+#         PARAMETER_ID: parameter_dict0['id'],
+#         PARAMETER_NAME: parameter_dict0['name'],
+#         PARAMETER_SCALE: LIN,
+#         NOMINAL_VALUE: 0,
+#         ESTIMATE: 0,
+#     }
+#     parameter_dicts.append(parameter_dict)
+## Basal
+for parameter_dict0 in common_parameters.values():
+    if parameter_dict0['id'] == 'a_b':
+        parameter_dict = {
+            PARAMETER_ID: parameter_dict0['id'],
+            PARAMETER_NAME: parameter_dict0['name'],
+            PARAMETER_SCALE: LOG10,
+            LOWER_BOUND: '1e-12',
+            UPPER_BOUND: '1e3',
+            NOMINAL_VALUE: 0,
+            ESTIMATE: 1,
+        }
+    elif parameter_dict0['id'] == 'da_b':
+        parameter_dict = {
+            PARAMETER_ID: parameter_dict0['id'],
+            PARAMETER_NAME: parameter_dict0['name'],
+            PARAMETER_SCALE: LIN,
+            NOMINAL_VALUE: 1,
+            ESTIMATE: 0,
+        }
+    else:
+        raise NotImplementedError(parameter_dict0['id'])
+    parameter_dicts.append(parameter_dict)
+## Noise
+parameter_dicts.append(
+    {
+        PARAMETER_ID: noise,
+        PARAMETER_NAME: noise,
+        PARAMETER_SCALE: LOG10,
+        LOWER_BOUND: '1e-12',
+        UPPER_BOUND: '1e3',
+        NOMINAL_VALUE: 0.1,
+        ESTIMATE: 1,
+    }
+)
+
+condition_df = petab.get_condition_df(pd.DataFrame({CONDITION_ID:[condition_id]}))
+observable_df = petab.get_observable_df(pd.DataFrame(observable_dicts))
+measurement_df = petab.get_measurement_df(pd.DataFrame(measurement_dicts))
+parameter_df = petab.get_parameter_df(pd.DataFrame(parameter_dicts))
+
+petab.write_condition_df(condition_df,     'output/petab/conditions.tsv')
+petab.write_observable_df(observable_df,   'output/petab/observables.tsv')
+petab.write_measurement_df(measurement_df, 'output/petab/measurements.tsv')
+petab.write_parameter_df(parameter_df,     'output/petab/parameters.tsv')
+shutil.copy('input/petab_problem.yaml',    'output/petab/petab_problem.yaml')

Benchmark-Models/Blasi_CellSystems2016/petab_select/04_create_select.py

@@ -0,0 +1,35 @@
+from model_info import (
+    acetylation_reactions,
+    get_parameter_id,
+)
+import shutil
+import pandas as pd
+
+
+import petab_select
+from petab_select.constants import (
+    MODEL_SUBSPACE_ID,
+    PETAB_YAML,
+    ESTIMATE,
+    PARAMETER_VALUE_DELIMITER,
+)
+
+
+parameter_ids = [get_parameter_id(*reaction) for reaction in acetylation_reactions]
+
+# The value to "turn off" the parameter. Since it's a multiplier, `1` is the identity.
+default_motif_value = '1.0'
+
+model_subspace = {
+    MODEL_SUBSPACE_ID: 'M',
+    PETAB_YAML: '../petab/petab_problem.yaml',
+    **{
+        parameter_id: PARAMETER_VALUE_DELIMITER.join([default_motif_value, ESTIMATE])
+        for parameter_id in parameter_ids
+    }
+}
+
+df = petab_select.get_model_space_df(pd.DataFrame([model_subspace]))
+
+petab_select.write_model_space_df(df, filename='output/select/model_space.tsv')
+shutil.copy('input/petab_select_problem.yaml', 'output/select/petab_select_problem.yaml')

Benchmark-Models/Blasi_CellSystems2016/petab_select/README.md

@@ -0,0 +1,11 @@
+# Model selection problem for `Blasi_CellSystems2016`
+This is a fresh implementation of the model from `Blasi_CellSystems2016`, to faciliate model selection with PEtab and PEtab Select. The implementation can be found in the `output/` directory.
+
+Specifically, this model contains parameters to enable model selection with motif-specific acetylation rates. These additional parameters are implemented as multipliers of the basal acetylation rate `a_b`.
+Note that this means any motif-specific parameter `a_*` (except the basal acetylation rate `a_b`) can be fixed to 1 to "turn off" the motif-specific acetylation rate, or estimated to "turn on" the motif-specific acetylation rate.
+
+The only thing taken from the main PEtab problem are the measurements, the remaining PEtab (Select) files and the model were generated with the provided scripts.
+
+# Python requirements
+The following should be sufficient to run the provided scripts. If using a system that doesn't support the bash `*.sh` script, then simply create the folders manually.
+- a Python 3 virtual environment with PEtab Select (`pip install petab-select`) installed.

Benchmark-Models/Blasi_CellSystems2016/petab_select/input/petab_problem.yaml

@@ -0,0 +1,11 @@
+format_version: 1
+parameter_file: parameters.tsv
+problems:
+- condition_files:
+  - conditions.tsv
+  measurement_files:
+  - measurements.tsv
+  observable_files:
+  - observables.tsv
+  sbml_files:
+  - ../model/model.xml

Benchmark-Models/Blasi_CellSystems2016/petab_select/input/petab_select_problem.yaml

@@ -0,0 +1,5 @@
+format_version: beta_1
+criterion: BIC
+method: forward
+model_space_files:
+- model_space.tsv