Merge pull request #422 from biosustain/feat-conc-init-exp

Feat: specify initial concentrations individually for the ODE state vector

docs/usage/inputting.rst

@@ -333,6 +333,7 @@ experiment.
 - :code:`pme_knockouts`: An optional table describing knockouts of
   phosphorylation modifying enzymes. Each entry has one field called "pme".
 - :code:`measurements`: An optional table describing measurements
+- :code:`initial_state`: An optional table describing the initial concentrations of balanced metabolites
 
 The measurement table has these fields:
 
@@ -353,6 +354,19 @@ The measurement table has these fields:
 measurements or the scale parameter of a lognormal distribution for
 concentration measurements.
 
+The :code:`initial_state` table has these fields:
+
+- :code:`metabolite` A string identifying the metabolite.
+- :code:`compartment` A string identifying the compartment.
+- :code:`value` The initial state, as a float.
+
+A default initial concentration for all metabolites and experiments can be
+specified in the general configuration. In addition, this :code:`initial_state`
+table can be specified for a more granular control of the initial vector of
+concentrations to solve the steady state ODE. This may be useful when the
+concentrations in the model span over several orders of magnitudes and thus a
+unique initial value is suboptimal.
+
 The priors file
 ===============
 

src/maud/data/example_inputs/example_ode/experiments.toml

@@ -7,4 +7,7 @@ measurements = [
   {target_type = "mic", metabolite = "C", compartment = "c", value = 3.021870, error_scale = 0.1},
   {target_type = "flux", reaction = "r3", value = 0.421816, error_scale = 0.01},
   {target_type = "flux", reaction = "r4", value = 2.116740, error_scale = 0.01},
-] 
+]
+initial_state = [
+  {metabolite = "D", compartment = "c", value = 3.0},
+]

src/maud/data/example_inputs/linear/expected_stan_input.json

@@ -1 +1 @@
-{"N_mic": 4, "N_edge_sub": 3, "N_edge_prod": 3, "N_edge": 3, "N_unbalanced": 2, "N_enzyme": 3, "N_phosphorylation": 0, "N_pme": 0, "N_competitive_inhibition": 1, "N_allostery": 2, "N_allosteric_enzyme": 2, "N_drain": 0, "N_km": 5, "N_sub_km": 3, "N_prod_km": 2, "S": [[-1, 0, 0], [1, -1, 0], [0, 1, -1], [0, 0, 1]], "N_reaction": 3, "N_metabolite": 2, "balanced_mic_ix": [2, 3], "unbalanced_mic_ix": [1, 4], "ci_mic_ix": [2], "edge_type": [1, 2, 1], "edge_to_enzyme": [1, 2, 3], "edge_to_tc": [1, 2, 0], "edge_to_drain": [0, 0, 0], "edge_to_reaction": [1, 2, 3], "water_stoichiometry": [0.0, 0.0, 0.0], "transported_charge": [0.0, 0.0, 1.0], "mic_to_met": [1, 1, 2, 2], "subunits": [1, 1, 1], "sub_by_edge_long": [1, 2, 3], "sub_by_edge_bounds": [[1, 1], [2, 2], [3, 3]], "prod_by_edge_long": [2, 3, 4], "prod_by_edge_bounds": [[1, 1], [2, 2], [3, 3]], "sub_km_ix_by_edge_long": [1, 3, 4], "sub_km_ix_by_edge_bounds": [[1, 1], [2, 2], [3, 3]], "prod_km_ix_by_edge_long": [2, 5], "prod_km_ix_by_edge_bounds": [[1, 1], [2, 1], [2, 2]], "ci_ix_long": [1], "ci_ix_bounds": [[1, 0], [1, 1], [2, 1]], "allostery_ix_long": [1, 2], "allostery_ix_bounds": [[1, 1], [2, 2], [3, 2]], "allostery_type": [1, 2], "allostery_mic": [3, 2], "phosphorylation_ix_long": [], "phosphorylation_ix_bounds": [[1, 0], [1, 0], [1, 0]], "phosphorylation_type": [], "phosphorylation_pme": [], "priors_km": [[0.0, 0.0, 0.0, 0.0, 0.0], [0.6, 0.6, 0.6, 0.6, 0.6]], "priors_ki": [[0.0], [0.6]], "priors_kcat": [[0.0, 0.0, 0.0], [0.6, 0.6, 0.6]], "priors_dissociation_constant": [[0.0, 0.0], [0.6, 0.6]], "priors_transfer_constant": [[0.0, 0.0], [0.6, 0.6]], "priors_kcat_pme": [[], []], "priors_drain_train": [[[]], [[]]], "priors_conc_enzyme_train": [[[0.5, 0.5, 0.5], [0.5, 0.5, 0.5]], [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]], "priors_conc_unbalanced_train": [[[-2.3, -2.3], [-2.3, -2.3]], [[2.0, 2.0], [2.0, 2.0]]], "priors_conc_pme_train": [[[]], [[]]], "priors_psi_train": [[-0.95, -0.95], [0.2, 0.2]], "prior_loc_dgf": [-1.0, -2.0], "prior_cov_dgf": [[0.05, 0.0], [0.0, 0.05]], "N_experiment_train": 2, "N_flux_measurement_train": 2, "N_enzyme_measurement_train": 6, "N_conc_measurement_train": 8, "N_enzyme_knockout_train": 0, "N_pme_knockout_train": 0, "temperature_train": [299.0, 298.15], "enzyme_knockout_train_long": [], "enzyme_knockout_train_bounds": [[1, 0], [1, 0]], "pme_knockout_train_long": [], "pme_knockout_train_bounds": [[1, 0], [1, 0]], "yconc_train": [0.59, 0.38, 1.09, 1.05, 0.54, 0.38, 1.12, 1.14], "sigma_yconc_train": [0.1, 0.1, 0.05, 0.05, 0.1, 0.1, 0.05, 0.05], "experiment_yconc_train": [1, 1, 1, 1, 2, 2, 2, 2], "mic_ix_yconc_train": [2, 3, 1, 4, 2, 3, 1, 4], "yflux_train": [0.19, 0.39], "sigma_yflux_train": [0.1, 0.1], "experiment_yflux_train": [1, 2], "reaction_yflux_train": [3, 3], "yenz_train": [1.5, 1.5, 1.5, 1.5, 1.5, 1.5], "sigma_yenz_train": [0.1, 0.1, 0.1, 0.1, 0.1, 0.1], "experiment_yenz_train": [1, 1, 1, 2, 2, 2], "enzyme_yenz_train": [1, 2, 3, 1, 2, 3], "likelihood": 1, "drain_small_conc_corrector": 1e-06, "reject_non_steady": 1, "steady_state_threshold_abs": 1e-06, "steady_state_threshold_rel": 0.001, "rel_tol": 1e-05, "abs_tol": 1e-05, "timepoint": 1000.0, "max_num_steps": 1000000, "conc_init": [[0.59, 0.38], [0.54, 0.38]]}
+{"N_mic": 4, "N_edge_sub": 3, "N_edge_prod": 3, "N_edge": 3, "N_unbalanced": 2, "N_enzyme": 3, "N_phosphorylation": 0, "N_pme": 0, "N_competitive_inhibition": 1, "N_allostery": 2, "N_allosteric_enzyme": 2, "N_drain": 0, "N_km": 5, "N_sub_km": 3, "N_prod_km": 2, "S": [[-1, 0, 0], [1, -1, 0], [0, 1, -1], [0, 0, 1]], "N_reaction": 3, "N_metabolite": 2, "balanced_mic_ix": [2, 3], "unbalanced_mic_ix": [1, 4], "ci_mic_ix": [2], "edge_type": [1, 2, 1], "edge_to_enzyme": [1, 2, 3], "edge_to_tc": [1, 2, 0], "edge_to_drain": [0, 0, 0], "edge_to_reaction": [1, 2, 3], "water_stoichiometry": [0.0, 0.0, 0.0], "transported_charge": [0.0, 0.0, 1.0], "mic_to_met": [1, 1, 2, 2], "subunits": [1, 1, 1], "sub_by_edge_long": [1, 2, 3], "sub_by_edge_bounds": [[1, 1], [2, 2], [3, 3]], "prod_by_edge_long": [2, 3, 4], "prod_by_edge_bounds": [[1, 1], [2, 2], [3, 3]], "sub_km_ix_by_edge_long": [1, 3, 4], "sub_km_ix_by_edge_bounds": [[1, 1], [2, 2], [3, 3]], "prod_km_ix_by_edge_long": [2, 5], "prod_km_ix_by_edge_bounds": [[1, 1], [2, 1], [2, 2]], "ci_ix_long": [1], "ci_ix_bounds": [[1, 0], [1, 1], [2, 1]], "allostery_ix_long": [1, 2], "allostery_ix_bounds": [[1, 1], [2, 2], [3, 2]], "allostery_type": [1, 2], "allostery_mic": [3, 2], "phosphorylation_ix_long": [], "phosphorylation_ix_bounds": [[1, 0], [1, 0], [1, 0]], "phosphorylation_type": [], "phosphorylation_pme": [], "priors_km": [[0.0, 0.0, 0.0, 0.0, 0.0], [0.6, 0.6, 0.6, 0.6, 0.6]], "priors_ki": [[0.0], [0.6]], "priors_kcat": [[0.0, 0.0, 0.0], [0.6, 0.6, 0.6]], "priors_dissociation_constant": [[0.0, 0.0], [0.6, 0.6]], "priors_transfer_constant": [[0.0, 0.0], [0.6, 0.6]], "priors_kcat_pme": [[], []], "priors_drain_train": [[[]], [[]]], "priors_conc_enzyme_train": [[[0.5, 0.5, 0.5], [0.5, 0.5, 0.5]], [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0]]], "priors_conc_unbalanced_train": [[[-2.3, -2.3], [-2.3, -2.3]], [[2.0, 2.0], [2.0, 2.0]]], "priors_conc_pme_train": [[[]], [[]]], "priors_psi_train": [[-0.95, -0.95], [0.2, 0.2]], "prior_loc_dgf": [-1.0, -2.0], "prior_cov_dgf": [[0.05, 0.0], [0.0, 0.05]], "N_experiment_train": 2, "N_flux_measurement_train": 2, "N_enzyme_measurement_train": 6, "N_conc_measurement_train": 7, "N_enzyme_knockout_train": 0, "N_pme_knockout_train": 0, "temperature_train": [299.0, 298.15], "enzyme_knockout_train_long": [], "enzyme_knockout_train_bounds": [[1, 0], [1, 0]], "pme_knockout_train_long": [], "pme_knockout_train_bounds": [[1, 0], [1, 0]], "yconc_train": [0.59, 1.09, 1.05, 0.54, 0.38, 1.12, 1.14], "sigma_yconc_train": [0.1, 0.05, 0.05, 0.1, 0.1, 0.05, 0.05], "experiment_yconc_train": [1, 1, 1, 2, 2, 2, 2], "mic_ix_yconc_train": [2, 1, 4, 2, 3, 1, 4], "yflux_train": [0.19, 0.39], "sigma_yflux_train": [0.1, 0.1], "experiment_yflux_train": [1, 2], "reaction_yflux_train": [3, 3], "yenz_train": [1.5, 1.5, 1.5, 1.5, 1.5, 1.5], "sigma_yenz_train": [0.1, 0.1, 0.1, 0.1, 0.1, 0.1], "experiment_yenz_train": [1, 1, 1, 2, 2, 2], "enzyme_yenz_train": [1, 2, 3, 1, 2, 3], "likelihood": 1, "drain_small_conc_corrector": 1e-06, "reject_non_steady": 1, "steady_state_threshold_abs": 1e-06, "steady_state_threshold_rel": 0.001, "rel_tol": 1e-05, "abs_tol": 1e-05, "timepoint": 1000.0, "max_num_steps": 1000000, "conc_init": [[0.59, 0.38], [0.54, 0.38]]}

src/maud/data/example_inputs/linear/experiments.toml

@@ -5,14 +5,16 @@ is_test = true
 temperature = 299.0
 measurements = [
   {target_type = "mic", metabolite = "M1", compartment = "c", value = 0.59, error_scale = 0.1},
-  {target_type = "mic", metabolite = "M2", compartment = "c", value = 0.38, error_scale = 0.1},
   {target_type = "mic", metabolite = "M1", compartment = "e", value = 1.09, error_scale = 0.05},
   {target_type = "mic", metabolite = "M2", compartment = "e", value = 1.05, error_scale = 0.05},
   {target_type = "flux", reaction = "r3", value = 0.19, error_scale = 0.1},
   {target_type = "enzyme", enzyme = "r1", value = 1.5, error_scale = 0.1},
   {target_type = "enzyme", enzyme = "r2", value = 1.5, error_scale = 0.1},
   {target_type = "enzyme", enzyme = "r3", value = 1.5, error_scale = 0.1},
 ]
+initial_state = [
+  {metabolite = "M2", compartment = "c", value = 0.38},
+]
 
 [[experiment]]
 id = "condition2"

src/maud/data_model/experiment.py

@@ -75,6 +75,20 @@ def __post_init__(self):
         self.id = ID_SEPARATOR.join(["pko", self.experiment, self.enzyme])
 
 
+@dataclass
+class InitConcentration:
+    """Indication of the initial value of a concentration in the ODE."""
+
+    metabolite: str
+    compartment: str
+    value: float
+    target_id: str = Field(default=None, init=False, exclude=True)
+
+    def __post_init__(self):
+        """Add target_id field."""
+        self.target_id = ID_SEPARATOR.join([self.metabolite, self.compartment])
+
+
 @dataclass
 class Experiment:
     """Maud representation of an experiment.
@@ -89,6 +103,7 @@ class Experiment:
     is_test: bool
     temperature: float = 298.15
     measurements: List[Measurement] = Field(default_factory=lambda: [])
+    initial_state: List[InitConcentration] = Field(default_factory=lambda: [])
     enzyme_knockouts: List[EnzymeKnockout] = Field(default_factory=lambda: [])
     pme_knockouts: List[PhosphorylationModifyingEnzymeKnockout] = Field(
         default_factory=lambda: []

src/maud/getting_stan_inputs.py

@@ -509,7 +509,7 @@ def get_conc_init(
             m
             for m in experiment.measurements
             if m.target_type == MeasurementType.MIC
-        ]
+        ] + experiment.initial_state
         for j, mic in enumerate(balanced_mics):
             if any(mic.id == m.target_id for m in msmts):
                 inits[i][j] = gmean(