From 1a749d399a3931d94d98738a1b3b451cf4defb64 Mon Sep 17 00:00:00 2001 From: NicholasCowie Date: Mon, 24 Jun 2024 14:32:44 +0200 Subject: [PATCH] feat: onserved moieties (#459) * feat: left_nullspace for conserved moieties * feat: including data structure and assertions for conserved moieties * feat: data structure for conserved moieties * change np array to df * feat: initialised at fraction of sampled pool size and added ODE solver * bug: error in order of conserved moiety groups * fix: add change penalty term * fix: change penalty term to abs based on #434 * feat: include example of conserved moieties * feat: including printout of moeity pools and sorting check * fixing out_of_sample_model * including corrected arviz branch, fix next arviz release * including corrected arviz branch, fix next arviz release * fix: sort properly the conserved_moieties to match their existence * bug: correcting ordering of conserved moiety pool indexes * test: fixing example odes and moving new ode example * test: include json files * chore: cleanup new conserved moiety example --------- Co-authored-by: carrascomj --- maud/cli.py | 9 + .../example_ode/experimental_setup.toml | 4 - .../example_inputs/example_ode/inits.json | 2 +- .../example_ode/input_data_train.json | 2 +- .../{ => example_ode_2}/__init__.py | 0 .../example_inputs/example_ode_2/config.toml | 22 + .../example_ode_2/experiments.toml | 189 ++++++ .../example_ode_2/kinetic_model.toml | 159 ++++++ .../example_inputs/example_ode_2/priors.toml | 205 +++++++ .../linear/expected_stan_input.json | 536 +----------------- maud/data_model/experiment.py | 9 +- maud/data_model/id_component.py | 1 + maud/data_model/kinetic_model.py | 76 +++ maud/data_model/maud_config.py | 5 +- maud/data_model/maud_init.py | 1 + maud/data_model/maud_parameter.py | 14 + maud/data_model/parameter_input.py | 4 + maud/data_model/parameter_set.py | 26 + maud/getting_idatas.py | 21 +- maud/getting_stan_inputs.py | 66 ++- maud/parsing_kinetic_models.py | 13 + maud/running_stan.py | 4 +- maud/stan/functions.stan | 59 +- maud/stan/model.stan | 134 +++-- maud/stan/out_of_sample_model.stan | 93 +-- maud/utility_functions.py | 48 ++ pyproject.toml | 10 +- 27 files changed, 1028 insertions(+), 684 deletions(-) delete mode 100644 maud/data/example_inputs/example_ode/experimental_setup.toml rename maud/data/example_inputs/{ => example_ode_2}/__init__.py (100%) create mode 100644 maud/data/example_inputs/example_ode_2/config.toml create mode 100644 maud/data/example_inputs/example_ode_2/experiments.toml create mode 100644 maud/data/example_inputs/example_ode_2/kinetic_model.toml create mode 100644 maud/data/example_inputs/example_ode_2/priors.toml diff --git a/maud/cli.py b/maud/cli.py index 216062e61..87781f2fa 100644 --- a/maud/cli.py +++ b/maud/cli.py @@ -219,6 +219,15 @@ def do_simulate(data_path, output_dir, n): .unstack() .T ) + print("\n\nSimulated conserved moiety pools:") + if "conc_moiety_pool_train" in idata.posterior.keys(): + print( + idata.posterior["conc_moiety_pool_train"] + .mean(dim=["chain", "draw"]) + .to_series() + .unstack() + .T + ) print("\n\nSimulated fluxes:") print( idata.posterior["flux_train"] diff --git a/maud/data/example_inputs/example_ode/experimental_setup.toml b/maud/data/example_inputs/example_ode/experimental_setup.toml deleted file mode 100644 index fe8289b08..000000000 --- a/maud/data/example_inputs/example_ode/experimental_setup.toml +++ /dev/null @@ -1,4 +0,0 @@ -[[experiment]] -id = "condition1" -is_train = true -is_test = true diff --git a/maud/data/example_inputs/example_ode/inits.json b/maud/data/example_inputs/example_ode/inits.json index 96656dbe0..52c357fe8 100644 --- a/maud/data/example_inputs/example_ode/inits.json +++ b/maud/data/example_inputs/example_ode/inits.json @@ -1 +1 @@ -{"dgf": [10.0, 2.0, 5.0, 0.0], "dgf_free": [10.0, 2.0, 5.0, 0.0], "km": [0.5, 1.0, 2.0, 0.5, 1.0, 2.0, 0.5, 2.0, 1.0, 3.0000000000000004], "log_km_z": [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "ki": [1.0], "log_ki_z": [0.0], "kcat": [0.5, 2.0, 1.0, 2.0, 1.0], "log_kcat_z": [0.0, 0.0, 0.0, 0.0, 0.0], "dissociation_constant": [0.3, 0.9], "log_dissociation_constant_z": [0.0, 0.0], "transfer_constant": [1.0, 1.0], "log_transfer_constant_z": [0.0, 0.0], "kcat_pme": [], "log_kcat_pme_z": [], "drain_train": [[]], "drain_train_z": [[]], "drain_test": [[]], "drain_test_z": [[]], "conc_enzyme_train": [[1.0, 2.0, 3.0000000000000004, 2.0, 3.0000000000000004]], "log_conc_enzyme_train_z": [[0.0, 0.0, 0.0, 0.0, 0.0]], "conc_enzyme_test": [[1.0, 2.0, 3.0000000000000004, 2.0, 3.0000000000000004]], "log_conc_enzyme_test_z": [[0.0, 0.0, 0.0, 0.0, 0.0]], "conc_unbalanced_train": [[4.999999999999999, 0.5]], "log_conc_unbalanced_train_z": [[0.0, 0.0]], "conc_unbalanced_test": [[4.999999999999999, 0.5]], "log_conc_unbalanced_test_z": [[0.0, 0.0]], "conc_pme_train": [[]], "log_conc_pme_train_z": [[]], "conc_pme_test": [[]], "log_conc_pme_test_z": [[]], "psi_train": [0.0], "psi_train_z": [0.0], "psi_test": [0.0], "psi_test_z": [0.0]} +{"dgf": [10.0, 2.0, 5.0, 0.0], "dgf_free": [10.0, 2.0, 5.0, 0.0], "km": [0.5, 1.0, 2.0, 0.5, 1.0, 2.0, 0.5, 2.0, 1.0, 3.0000000000000004], "log_km_z": [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "ki": [1.0], "log_ki_z": [0.0], "kcat": [0.5, 2.0, 1.0, 2.0, 1.0], "log_kcat_z": [0.0, 0.0, 0.0, 0.0, 0.0], "dissociation_constant": [0.3, 0.9], "log_dissociation_constant_z": [0.0, 0.0], "transfer_constant": [1.0, 1.0], "log_transfer_constant_z": [0.0, 0.0], "kcat_pme": [], "log_kcat_pme_z": [], "drain_train": [[]], "drain_train_z": [[]], "drain_test": [[]], "drain_test_z": [[]], "conc_enzyme_train": [[1.0, 2.0, 3.0000000000000004, 2.0, 3.0000000000000004]], "log_conc_enzyme_train_z": [[0.0, 0.0, 0.0, 0.0, 0.0]], "conc_enzyme_test": [[1.0, 2.0, 3.0000000000000004, 2.0, 3.0000000000000004]], "log_conc_enzyme_test_z": [[0.0, 0.0, 0.0, 0.0, 0.0]], "conc_unbalanced_train": [[4.999999999999999, 0.5]], "log_conc_unbalanced_train_z": [[0.0, 0.0]], "conc_unbalanced_test": [[4.999999999999999, 0.5]], "log_conc_unbalanced_test_z": [[0.0, 0.0]], "conc_pme_train": [[]], "log_conc_pme_train_z": [[]], "conc_pme_test": [[]], "log_conc_pme_test_z": [[]], "psi_train": [0.0], "psi_train_z": [0.0], "psi_test": [0.0], "psi_test_z": [0.0]} \ No newline at end of file diff --git a/maud/data/example_inputs/example_ode/input_data_train.json b/maud/data/example_inputs/example_ode/input_data_train.json index 40a4a01ad..48ecacac1 100644 --- a/maud/data/example_inputs/example_ode/input_data_train.json +++ b/maud/data/example_inputs/example_ode/input_data_train.json @@ -1 +1 @@ -{"N_mic": 4, "N_edge_sub": 5, "N_edge_prod": 5, "N_edge": 5, "N_unbalanced": 2, "N_enzyme": 5, "N_er": 5, "N_phosphorylation": 0, "N_pme": 0, "N_competitive_inhibition": 1, "N_allostery": 2, "N_allosteric_enzyme": 2, "N_drain": 0, "N_km": 10, "N_sub_km": 5, "N_prod_km": 5, "S": [[-1, -1, -1, 0, 0], [1, 0, 0, -1, 0], [0, 1, 1, 0, -1], [0, 0, 0, 1, 1]], "N_reaction": 4, "N_metabolite": 4, "balanced_mic_ix": [2, 3], "unbalanced_mic_ix": [1, 4], "ci_mic_ix": [4], "edge_type": [1, 1, 1, 1, 1], "edge_to_enzyme": [1, 2, 3, 4, 5], "edge_to_er": [1, 2, 3, 4, 5], "edge_to_tc": [0, 1, 2, 0, 0], "edge_to_drain": [0, 0, 0, 0, 0], "edge_to_reaction": [1, 2, 2, 3, 4], "water_stoichiometry": [0.0, 0.0, 0.0, 0.0, 0.0], "transported_charge": [0.0, 0.0, 0.0, 0.0, 0.0], "mic_to_met": [1, 2, 3, 4], "subunits": [1, 1, 1, 1, 1], "sub_by_edge_long": [1, 1, 1, 2, 3], "sub_by_edge_bounds": [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]], "prod_by_edge_long": [2, 3, 3, 4, 4], "prod_by_edge_bounds": [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]], "sub_km_ix_by_edge_long": [1, 3, 5, 7, 9], "sub_km_ix_by_edge_bounds": [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]], "prod_km_ix_by_edge_long": [2, 4, 6, 8, 10], "prod_km_ix_by_edge_bounds": [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]], "ci_ix_long": [1], "ci_ix_bounds": [[1, 1], [2, 1], [2, 1], [2, 1], [2, 1]], "allostery_ix_long": [1, 2], "allostery_ix_bounds": [[1, 0], [1, 1], [2, 2], [3, 2], [3, 2]], "allostery_type": [1, 2], "allostery_mic": [3, 3], "phosphorylation_ix_long": [], "phosphorylation_ix_bounds": [[1, 0], [1, 0], [1, 0], [1, 0], [1, 0]], "phosphorylation_type": [], "phosphorylation_pme": [], "N_dgf_fixed": 0, "dgf_fixed": [], "ix_dgf_free": [1, 2, 3, 4], "ix_dgf_fixed": [], "priors_km": [[-0.6931471805599453, 0.0, 0.6931471805599453, -0.6931471805599453, 0.0, 0.6931471805599453, -0.6931471805599453, 0.6931471805599453, 0.0, 1.0986122886681098], [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2]], "priors_ki": [[0.0], [0.2]], "priors_kcat": [[-0.6931471805599453, 0.6931471805599453, 0.0, 0.6931471805599453, 0.0], [0.2, 0.2, 0.2, 0.2, 0.2]], "priors_dissociation_constant": [[-1.2039728043259361, -0.10536051565782628], [0.2, 0.2]], "priors_transfer_constant": [[0.0, 0.0], [0.2, 0.2]], "priors_kcat_pme": [[], []], "priors_drain_train": [[[]], [[]]], "priors_conc_enzyme_train": [[[0.0, 0.6931471805599453, 1.0986122886681098, 0.6931471805599453, 1.0986122886681098]], [[0.2, 0.2, 0.2, 0.2, 0.2]]], "priors_conc_unbalanced_train": [[[1.6094379124341003, -0.6931471805599453]], [[0.2, 0.2]]], "priors_conc_pme_train": [[[]], [[]]], "priors_psi_train": [[0.0], [2.0]], "prior_loc_dgf": [10.0, 2.0, 5.0, 0.0], "prior_cov_dgf": [[0.2, 0.0, 0.0, 0.0], [0.0, 0.2, 0.0, 0.0], [0.0, 0.0, 0.2, 0.0], [0.0, 0.0, 0.0, 0.2]], "N_experiment_train": 1, "N_flux_measurement_train": 2, "N_enzyme_measurement_train": 0, "N_conc_measurement_train": 2, "N_enzyme_knockout_train": 0, "N_pme_knockout_train": 0, "temperature_train": [298.15], "enzyme_knockout_train_long": [], "enzyme_knockout_train_bounds": [[1, 0]], "pme_knockout_train_long": [], "pme_knockout_train_bounds": [[1, 0]], "yconc_train": [0.323117, 3.02187], "sigma_yconc_train": [0.1, 0.1], "experiment_yconc_train": [1, 1], "mic_ix_yconc_train": [2, 3], "yflux_train": [0.421816, 2.11674], "sigma_yflux_train": [0.01, 0.01], "experiment_yflux_train": [1, 1], "reaction_yflux_train": [3, 4], "yenz_train": [], "sigma_yenz_train": [], "experiment_yenz_train": [], "enzyme_yenz_train": [], "likelihood": 1, "drain_small_conc_corrector": 1e-06, "reject_non_steady": 1, "penalize_non_steady": 0, "steady_state_threshold_abs": 1e-08, "steady_state_threshold_rel": 0.001, "steady_state_penalty_rel": 1e-08, "rel_tol_ode": 1e-09, "abs_tol_ode": 1e-09, "max_num_steps_ode": 1000000000, "rel_tol_alg": 1e-07, "abs_tol_alg": 1e-07, "max_num_steps_alg": 1000000, "conc_init": [[0.323117, 3.02187]]} +{"N_mic": 4, "N_pool": 0, "N_edge_sub": 5, "N_edge_prod": 5, "N_edge": 5, "N_unbalanced": 2, "N_independent": 2, "N_dependent": 0, "N_enzyme": 5, "N_er": 5, "N_phosphorylation": 0, "N_pme": 0, "N_competitive_inhibition": 1, "N_allostery": 2, "N_allosteric_enzyme": 2, "N_drain": 0, "N_km": 10, "N_sub_km": 5, "N_prod_km": 5, "S": [[-1, -1, -1, 0, 0], [1, 0, 0, -1, 0], [0, 1, 1, 0, -1], [0, 0, 0, 1, 1]], "left_nullspace_independent": [[]], "N_reaction": 4, "N_metabolite": 4, "independent_bal_ix": [2, 3], "dependent_bal_ix": [], "unbalanced_mic_ix": [1, 4], "ci_mic_ix": [4], "edge_type": [1, 1, 1, 1, 1], "edge_to_enzyme": [1, 2, 3, 4, 5], "edge_to_er": [1, 2, 3, 4, 5], "edge_to_tc": [0, 1, 2, 0, 0], "edge_to_drain": [0, 0, 0, 0, 0], "edge_to_reaction": [1, 2, 2, 3, 4], "water_stoichiometry": [0.0, 0.0, 0.0, 0.0, 0.0], "transported_charge": [0.0, 0.0, 0.0, 0.0, 0.0], "mic_to_met": [1, 2, 3, 4], "subunits": [1, 1, 1, 1, 1], "sub_by_edge_long": [1, 1, 1, 2, 3], "sub_by_edge_bounds": [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]], "prod_by_edge_long": [2, 3, 3, 4, 4], "prod_by_edge_bounds": [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]], "sub_km_ix_by_edge_long": [1, 3, 5, 7, 9], "sub_km_ix_by_edge_bounds": [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]], "prod_km_ix_by_edge_long": [2, 4, 6, 8, 10], "prod_km_ix_by_edge_bounds": [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]], "ci_ix_long": [1], "ci_ix_bounds": [[1, 1], [2, 1], [2, 1], [2, 1], [2, 1]], "allostery_ix_long": [1, 2], "allostery_ix_bounds": [[1, 0], [1, 1], [2, 2], [3, 2], [3, 2]], "allostery_type": [1, 2], "allostery_mic": [3, 3], "phosphorylation_ix_long": [], "phosphorylation_ix_bounds": [[1, 0], [1, 0], [1, 0], [1, 0], [1, 0]], "phosphorylation_type": [], "phosphorylation_pme": [], "N_dgf_fixed": 0, "dgf_fixed": [], "ix_dgf_free": [1, 2, 3, 4], "ix_dgf_fixed": [], "priors_km": [[-0.6931471805599453, 0.0, 0.6931471805599453, -0.6931471805599453, 0.0, 0.6931471805599453, -0.6931471805599453, 0.6931471805599453, 0.0, 1.0986122886681098], [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2]], "priors_ki": [[0.0], [0.2]], "priors_kcat": [[-0.6931471805599453, 0.6931471805599453, 0.0, 0.6931471805599453, 0.0], [0.2, 0.2, 0.2, 0.2, 0.2]], "priors_dissociation_constant": [[-1.2039728043259361, -0.10536051565782628], [0.2, 0.2]], "priors_transfer_constant": [[0.0, 0.0], [0.2, 0.2]], "priors_kcat_pme": [[], []], "priors_drain_train": [[[]], [[]]], "priors_conc_enzyme_train": [[[0.0, 0.6931471805599453, 1.0986122886681098, 0.6931471805599453, 1.0986122886681098]], [[0.2, 0.2, 0.2, 0.2, 0.2]]], "priors_conc_unbalanced_train": [[[1.6094379124341003, -0.6931471805599453]], [[0.2, 0.2]]], "priors_conc_pme_train": [[[]], [[]]], "priors_psi_train": [[0.0], [2.0]], "prior_loc_dgf": [10.0, 2.0, 5.0, 0.0], "prior_cov_dgf": [[0.2, 0.0, 0.0, 0.0], [0.0, 0.2, 0.0, 0.0], [0.0, 0.0, 0.2, 0.0], [0.0, 0.0, 0.0, 0.2]], "N_experiment_train": 1, "N_flux_measurement_train": 2, "N_enzyme_measurement_train": 0, "N_conc_measurement_train": 2, "N_enzyme_knockout_train": 0, "N_pme_knockout_train": 0, "temperature_train": [298.15], "enzyme_knockout_train_long": [], "enzyme_knockout_train_bounds": [[1, 0]], "pme_knockout_train_long": [], "pme_knockout_train_bounds": [[1, 0]], "yconc_train": [0.323117, 3.02187], "sigma_yconc_train": [0.1, 0.1], "experiment_yconc_train": [1, 1], "mic_ix_yconc_train": [2, 3], "yflux_train": [0.421816, 2.11674], "sigma_yflux_train": [0.01, 0.01], "experiment_yflux_train": [1, 1], "reaction_yflux_train": [3, 4], "yenz_train": [], "sigma_yenz_train": [], "experiment_yenz_train": [], "enzyme_yenz_train": [], "likelihood": 1, "drain_small_conc_corrector": 1e-06, "penalize_non_steady": 0, "steady_state_threshold_abs": 1e-08, "steady_state_threshold_rel": 0.001, "steady_state_penalty_abs": 1e-12, "rel_tol_ode": 1e-09, "abs_tol_ode": 1e-09, "timepoint": 10000.0, "max_num_steps_ode": 1000000000, "rel_tol_alg": 1e-07, "abs_tol_alg": 1e-07, "max_num_steps_alg": 1000000, "conc_init": [[0.323117, 3.02187]]} \ No newline at end of file diff --git a/maud/data/example_inputs/__init__.py b/maud/data/example_inputs/example_ode_2/__init__.py similarity index 100% rename from maud/data/example_inputs/__init__.py rename to maud/data/example_inputs/example_ode_2/__init__.py diff --git a/maud/data/example_inputs/example_ode_2/config.toml b/maud/data/example_inputs/example_ode_2/config.toml new file mode 100644 index 000000000..c9af5040e --- /dev/null +++ b/maud/data/example_inputs/example_ode_2/config.toml @@ -0,0 +1,22 @@ +name = "test-ode" +kinetic_model_file = "kinetic_model.toml" +priors_file = "priors.toml" +experiments_file = "experiments.toml" +likelihood = true + +[cmdstanpy_config] +iter_warmup = 200 +iter_sampling = 100 +chains = 4 +refresh = 1 +adapt_init_phase = 20 +metric = "dense_e" +adapt_metric_window = 5 + +[ode_config] +rel_tol = 1e-12 +abs_tol = 1e-12 +timepoint = 1e4 + +[quench_config] +quench_timepoints = [1e-8, 1e-6, 1e-4, 1e-2, 1e0, 1e1] diff --git a/maud/data/example_inputs/example_ode_2/experiments.toml b/maud/data/example_inputs/example_ode_2/experiments.toml new file mode 100644 index 000000000..29711fe20 --- /dev/null +++ b/maud/data/example_inputs/example_ode_2/experiments.toml @@ -0,0 +1,189 @@ +[[experiment]] +id = "condition1" +is_train = true +is_test = false +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.447, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 3.76, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.417, error_scale = 0.042 }, + { target_type = "flux", reaction = "r4", value = 0.750, error_scale = 0.075 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 0.5 }, + { metabolite = "C", compartment = "c", value = 3.5 }, +] + +[[experiment]] +id = "condition2" +is_train = true +is_test = false +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.271, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 2.44, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.423, error_scale = 0.042 }, + { target_type = "flux", reaction = "r4", value = 1.31, error_scale = 0.13 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 0.3 }, + { metabolite = "C", compartment = "c", value = 2.5 }, +] + +[[experiment]] +id = "condition3" +is_train = true +is_test = false +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.798, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 4.69, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.161, error_scale = 0.016 }, + { target_type = "flux", reaction = "r4", value = 0.553, error_scale = 0.055 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 0.8 }, + { metabolite = "C", compartment = "c", value = 4.7 }, +] + +[[experiment]] +id = "condition4" +is_train = true +is_test = false +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.461, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 3.50, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.875, error_scale = 0.088 }, + { target_type = "flux", reaction = "r4", value = 0.736, error_scale = 0.074 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 0.5 }, + { metabolite = "C", compartment = "c", value = 3.5 }, +] + +[[experiment]] +id = "condition5" +is_train = true +is_test = false +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.252, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 12.0, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.467, error_scale = 0.047 }, + { target_type = "flux", reaction = "r4", value = 0.453, error_scale = 0.045 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 0.3 }, + { metabolite = "C", compartment = "c", value = 12 }, +] + +[[experiment]] +id = "condition6" +is_train = false +is_test = true +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.361, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 11.4, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.242, error_scale = 0.024 }, + { target_type = "flux", reaction = "r4", value = 0.721, error_scale = 0.072 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 0.4 }, + { metabolite = "C", compartment = "c", value = 11 }, +] + +[[experiment]] +id = "condition7" +is_train = false +is_test = true +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.323, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 3.02, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.421, error_scale = 0.042 }, + { target_type = "flux", reaction = "r4", value = 2.11, error_scale = 0.21 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 0.3 }, + { metabolite = "C", compartment = "c", value = 3 }, +] + +[[experiment]] +id = "condition8" +is_train = false +is_test = true +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 7.98, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 1.35, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.721, error_scale = 0.072 }, + { target_type = "flux", reaction = "r4", value = 2.04, error_scale = 0.20 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 8 }, + { metabolite = "C", compartment = "c", value = 1 }, +] + +[[experiment]] +id = "condition9" +is_train = false +is_test = true +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.377, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 7.8, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.462, error_scale = 0.046 }, + { target_type = "flux", reaction = "r4", value = 0.690, error_scale = 0.069 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 4 }, + { metabolite = "C", compartment = "c", value = 8 }, +] + +[[experiment]] +id = "condition10" +is_train = false +is_test = true +temperature = 298 + +measurements = [ + { target_type = "mic", metabolite = "B", compartment = "c", value = 0.447, error_scale = 0.1 }, + { target_type = "mic", metabolite = "C", compartment = "c", value = 3.76, error_scale = 0.1 }, + { target_type = "flux", reaction = "r3", value = 0.209, error_scale = 0.021 }, + { target_type = "flux", reaction = "r4", value = 0.375, error_scale = 0.038 }, + { target_type = "conserved_moiety", moiety_group = ["X1_c", "X2_c"], moiety_id = "X", value = 1, error_scale = 0.1} +] + +initial_state = [ + { metabolite = "B", compartment = "c", value = 0.45 }, + { metabolite = "C", compartment = "c", value = 4 }, +] diff --git a/maud/data/example_inputs/example_ode_2/kinetic_model.toml b/maud/data/example_inputs/example_ode_2/kinetic_model.toml new file mode 100644 index 000000000..b64364259 --- /dev/null +++ b/maud/data/example_inputs/example_ode_2/kinetic_model.toml @@ -0,0 +1,159 @@ +compartment = [ + {id = "c", name = "cytosol", volume = 1} + {id = "e", name = "external", volume = 1} +] + +metabolite_in_compartment = [ + {metabolite_id = "A", compartment_id = "c", balanced = true}, + {metabolite_id = "A", compartment_id = "e", balanced = false}, + {metabolite_id = "B", compartment_id = "c", balanced = true}, + {metabolite_id = "C", compartment_id = "c", balanced = true}, + {metabolite_id = "D", compartment_id = "c", balanced = true}, + {metabolite_id = "D", compartment_id = "e", balanced = false}, + {metabolite_id = "X1", compartment_id = "c", balanced = true}, + {metabolite_id = "X2", compartment_id = "c", balanced = true}, + {metabolite_id = "Z", compartment_id = "c", balanced = false}, +] + +enzyme_reaction = [ + {enzyme_id = "transA", reaction_id = "transA"}, + {enzyme_id = "transD", reaction_id = "transD"}, + {enzyme_id = "r1", reaction_id = "r1"}, + {enzyme_id = "r2A", reaction_id = "r2"}, + {enzyme_id = "r2B", reaction_id = "r2"}, + {enzyme_id = "r3", reaction_id = "r3"}, + {enzyme_id = "r4", reaction_id = "r4"}, + {enzyme_id = "regX", reaction_id = "regX"}, +] + +[[conserved_moiety]] +id = "X" +name = "Conserved X" +moiety_group = ["X1_c", "X2_c"] + +[[metabolite]] +id = "A" +name = "Metabolite A" +inchi_key = "" + +[[metabolite]] +id = "B" +name = "Metabolite B" +inchi_key = "" + +[[metabolite]] +id = "C" +name = "Metabolite C" +inchi_key = "" + +[[metabolite]] +id = "D" +name = "Metabolite D" +inchi_key = "" + +[[metabolite]] +id = "X1" +name = "Cofactor X1" +inchi_key = "" + +[[metabolite]] +id = "X2" +name = "Cofactor X2" +inchi_key = "" + +[[metabolite]] +id = "Z" +name = "Cofactor Z" +inchi_key = "" + +[[enzyme]] +id = "transA" +name = 'Exchanges A metabolite between internal and external' + +[[enzyme]] +id = "transD" +name = 'Exchanges D metabolite between internal and external' + +[[enzyme]] +id = "regX" +name = 'Regenerates the X cofactor' + +[[enzyme]] +id = "r1" +name = 'the enzyme that catalyses reaction r1' + +[[enzyme]] +id = 'r2A' +name = 'Isoenzyme A that catalyses reaction r2' + +[[enzyme]] +id = 'r2B' +name = 'Isoenzyme B that catalyses reaction r2' + +[[enzyme]] +id = 'r3' +name = 'the enzyme that catalyses reaction r3' + +[[enzyme]] +id = 'r4' +name = 'the enzyme that catalyses reaction r4' + +[[reaction]] +id = "transA" +name = "transport of A" +stoichiometry = {A_e = -1, A_c = 1} +mechanism = "reversible_michaelis_menten" + +[[reaction]] +id = "transD" +name = "transport of D" +stoichiometry = {D_c = -1, D_e = 1} +mechanism = "reversible_michaelis_menten" + +[[reaction]] +id = "regX" +name = "regenerates X" +stoichiometry = {Z_c = -1, X2_c = -1, X1_c = 1} +mechanism = "reversible_michaelis_menten" + +[[reaction]] +id = "r1" +name = "reaction number 1" +stoichiometry = {A_c = -1, B_c = 1} +mechanism = "reversible_michaelis_menten" + +[[reaction]] +id = "r2" +name = "reaction number 2" +stoichiometry = {A_c = -1, C_c = 1} +mechanism = "reversible_michaelis_menten" + +[[reaction]] +id = "r3" +name = "reaction number 3" +stoichiometry = {X1_c = -1, B_c = -1, D_c = 1, X2_c = 1} +mechanism = "reversible_michaelis_menten" + +[[reaction]] +id = "r4" +name = "reaction number 4" +stoichiometry = {C_c = -1, D_c = 1} +mechanism = "reversible_michaelis_menten" + +[[competitive_inhibition]] +enzyme_id = "r1" +reaction_id = "r1" +metabolite_id = "D" +compartment_id = "c" + +[[allostery]] +modification_type = "activation" +enzyme_id = "r2A" +metabolite_id = "C" +compartment_id = "c" + +[[allostery]] +modification_type = "inhibition" +enzyme_id = "r2B" +metabolite_id = "C" +compartment_id = "c" diff --git a/maud/data/example_inputs/example_ode_2/priors.toml b/maud/data/example_inputs/example_ode_2/priors.toml new file mode 100644 index 000000000..3d419c813 --- /dev/null +++ b/maud/data/example_inputs/example_ode_2/priors.toml @@ -0,0 +1,205 @@ + +km = [ + {metabolite = "A", compartment = "e", enzyme = "transA", exploc = 1, scale = 1}, + {metabolite = "A", compartment = "c", enzyme = "transA", exploc = 1, scale = 1}, + {metabolite = "D", compartment = "e", enzyme = "transD", exploc = 1, scale = 1}, + {metabolite = "D", compartment = "c", enzyme = "transD", exploc = 1, scale = 1}, + {metabolite = "Z", compartment = "c", enzyme = "regX", exploc = 0.1, scale = 1}, + {metabolite = "X1", compartment = "c", enzyme = "regX", exploc = 1, scale = 1}, + {metabolite = "X2", compartment = "c", enzyme = "regX", exploc = 0.01, scale = 1}, + {metabolite = "A", compartment = "c", enzyme = "r1", exploc = 0.5, scale = 1}, + {metabolite = "B", compartment = "c", enzyme = "r1", exploc = 1, scale = 1}, + {metabolite = "A", compartment = "c", enzyme = "r2A", exploc = 2, scale = 1}, + {metabolite = "C", compartment = "c", enzyme = "r2A", exploc = 0.5, scale = 1}, + {metabolite = "A", compartment = "c", enzyme = "r2B", exploc = 1, scale = 1}, + {metabolite = "C", compartment = "c", enzyme = "r2B", exploc = 2, scale = 1}, + {metabolite = "B", compartment = "c", enzyme = "r3", exploc = 0.5, scale = 1}, + {metabolite = "X2", compartment = "c", enzyme = "r3", exploc = 0.5, scale = 1}, + {metabolite = "X1", compartment = "c", enzyme = "r3", exploc = 0.01, scale = 1}, + {metabolite = "D", compartment = "c", enzyme = "r3", exploc = 2, scale = 1}, + {metabolite = "C", compartment = "c", enzyme = "r4", exploc = 1, scale = 1}, + {metabolite = "D", compartment = "c", enzyme = "r4", exploc = 3, scale = 1}, +] + +ki = [ + {metabolite = "D", compartment = "c", enzyme = "r1", reaction = "r1", exploc = 1, scale = 1}, +] + +dissociation_constant = [ + {metabolite = "C", compartment = "c", enzyme = "r2B", exploc = 0.9, scale = 1, modification_type = "inhibition"}, + {metabolite = "C", compartment = "c", enzyme = "r2A", exploc = 0.3, scale = 1, modification_type = "activation"}, +] + +transfer_constant = [ + {enzyme = "r2A", exploc = 1, scale = 1}, + {enzyme = "r2B", exploc = 1, scale = 1}, +] + +kcat = [ + {enzyme = "transA", reaction = "transA", exploc = 0.5, scale = 1}, + {enzyme = "transD", reaction = "transD", exploc = 2, scale = 1}, + {enzyme = "regX", reaction = "regX", exploc = 4, scale = 1}, + {enzyme = "r1", reaction = "r1", exploc = 0.5, scale = 1}, + {enzyme = "r2A", reaction = "r2", exploc = 2, scale = 1}, + {enzyme = "r2B", reaction = "r2", exploc = 1, scale = 1}, + {enzyme = "r3", reaction = "r3", exploc = 2, scale = 1}, + {enzyme = "r4", reaction = "r4", exploc = 1, scale = 1}, +] + +conc_unbalanced = [ + {metabolite = "A", compartment = "e", experiment = "condition1", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition1", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition1", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition2", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition2", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition2", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition3", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition3", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition3", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition4", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition4", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition4", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition5", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition5", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition5", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition6", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition6", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition6", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition7", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition7", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition7", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition8", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition8", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition8", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition9", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition9", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition9", exploc = 5, scale = 0.1}, + + {metabolite = "A", compartment = "e", experiment = "condition10", exploc = 5, scale = 0.1}, + {metabolite = "D", compartment = "e", experiment = "condition10", exploc = 0.5, scale = 0.1}, + {metabolite = "Z", compartment = "e", experiment = "condition10", exploc = 5, scale = 0.1}, +] + +conc_enzyme = [ + + {enzyme = "transA", experiment = "eondition1", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "eondition1", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition1", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition1", exploc = 1, scale = 0.1}, + {enzyme = "r2A", experiment = "condition1", exploc = 1, scale = 0.1}, + {enzyme = "r2B", experiment = "condition1", exploc = 1, scale = 0.1}, + {enzyme = "r3", experiment = "condition1", exploc = 1, scale = 0.1}, + {enzyme = "r4", experiment = "condition1", exploc = 1, scale = 0.1}, + + {enzyme = "transA", experiment = "condition2", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition2", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition2", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition2", exploc = 1, scale = 0.1}, + {enzyme = "r2A", experiment = "condition2", exploc = 2, scale = 0.1}, + {enzyme = "r2B", experiment = "condition2", exploc = 0.5, scale = 0.1}, + {enzyme = "r3", experiment = "condition2", exploc = 4, scale = 0.1}, + {enzyme = "r4", experiment = "condition2", exploc = 2, scale = 0.1}, + + {enzyme = "transA", experiment = "condition3", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition3", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition3", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition3", exploc = 0.4, scale = 0.1}, + {enzyme = "r2A", experiment = "condition3", exploc = 0.3, scale = 0.1}, + {enzyme = "r2B", experiment = "condition3", exploc = 1.5, scale = 0.1}, + {enzyme = "r3", experiment = "condition3", exploc = 0.2, scale = 0.1}, + {enzyme = "r4", experiment = "condition3", exploc = 0.7, scale = 0.1}, + + {enzyme = "transA", experiment = "condition4", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition4", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition4", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition4", exploc = 2.1, scale = 0.1}, + {enzyme = "r2A", experiment = "condition4", exploc = 1.3, scale = 0.1}, + {enzyme = "r2B", experiment = "condition4", exploc = 0.5, scale = 0.1}, + {enzyme = "r3", experiment = "condition4", exploc = 2, scale = 0.1}, + {enzyme = "r4", experiment = "condition4", exploc = 1, scale = 0.1}, + + {enzyme = "transA", experiment = "condition5", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition5", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition5", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition5", exploc = 1.1, scale = 0.1}, + {enzyme = "r2A", experiment = "condition5", exploc = 2.4, scale = 0.1}, + {enzyme = "r2B", experiment = "condition5", exploc = 1.2, scale = 0.1}, + {enzyme = "r3", experiment = "condition5", exploc = 7, scale = 0.1}, + {enzyme = "r4", experiment = "condition5", exploc = 0.5, scale = 0.1}, + + {enzyme = "transA", experiment = "condition6", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition6", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition6", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition6", exploc = 0.6, scale = 0.1}, + {enzyme = "r2A", experiment = "condition6", exploc = 1.2, scale = 0.1}, + {enzyme = "r2B", experiment = "condition6", exploc = 4.1, scale = 0.1}, + {enzyme = "r3", experiment = "condition6", exploc = 0.9, scale = 0.1}, + {enzyme = "r4", experiment = "condition6", exploc = 0.8, scale = 0.1}, + + {enzyme = "transA", experiment = "condition7", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition7", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition7", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition7", exploc = 1, scale = 0.1}, + {enzyme = "r2A", experiment = "condition7", exploc = 2, scale = 0.1}, + {enzyme = "r2B", experiment = "condition7", exploc = 3, scale = 0.1}, + {enzyme = "r3", experiment = "condition7", exploc = 2, scale = 0.1}, + {enzyme = "r4", experiment = "condition7", exploc = 3, scale = 0.1}, + + {enzyme = "transA", experiment = "condition8", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition8", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition8", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition8", exploc = 3, scale = 0.1}, + {enzyme = "r2A", experiment = "condition8", exploc = 2, scale = 0.1}, + {enzyme = "r2B", experiment = "condition8", exploc = 1, scale = 0.1}, + {enzyme = "r3", experiment = "condition8", exploc = 0.4, scale = 0.1}, + {enzyme = "r4", experiment = "condition8", exploc = 4, scale = 0.1}, + + {enzyme = "transA", experiment = "condition9", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition9", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition9", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition9", exploc = 1.1, scale = 0.1}, + {enzyme = "r2A", experiment = "condition9", exploc = 0.4, scale = 0.1}, + {enzyme = "r2B", experiment = "condition9", exploc = 3, scale = 0.1}, + {enzyme = "r3", experiment = "condition9", exploc = 1.5, scale = 0.1}, + {enzyme = "r4", experiment = "condition9", exploc = 0.8, scale = 0.1}, + + {enzyme = "transA", experiment = "condition10", exploc = 1, scale = 0.1}, + {enzyme = "transD", experiment = "condition10", exploc = 1, scale = 0.1}, + {enzyme = "regX", experiment = "condition10", exploc = 1, scale = 0.1}, + {enzyme = "r1", experiment = "condition10", exploc = 0.5, scale = 0.1}, + {enzyme = "r2A", experiment = "condition10", exploc = 0.5, scale = 0.1}, + {enzyme = "r2B", experiment = "condition10", exploc = 0.5, scale = 0.1}, + {enzyme = "r3", experiment = "condition10", exploc = 0.5, scale = 0.1}, + {enzyme = "r4", experiment = "condition10", exploc = 0.5, scale = 0.1}, +] + +dgf = [ + {metabolite = "A", location = 10.0, scale = 0.2}, + {metabolite = "B", location = 2.0, scale = 0.2}, + {metabolite = "C", location = 5.0, scale = 0.2}, + {metabolite = "D", location = 0.0, scale = 0.2}, + {metabolite = "X1", location = 1.0, scale = 0.2}, + {metabolite = "X2", location = 0.0, scale = 0.2}, + {metabolite = "Z", location = 5.0, scale = 0.2}, +] + +conc_moiety_pool = [ + {conserved_moiety= "X", experiment = "condition1", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition2", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition3", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition4", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition5", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition6", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition7", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition8", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition9", exploc = 2, scale = 0.1}, + {conserved_moiety= "X", experiment = "condition10", exploc = 2, scale = 0.1}, +] diff --git a/maud/data/example_inputs/linear/expected_stan_input.json b/maud/data/example_inputs/linear/expected_stan_input.json index 58b9105ba..b41b13423 100644 --- a/maud/data/example_inputs/linear/expected_stan_input.json +++ b/maud/data/example_inputs/linear/expected_stan_input.json @@ -1,535 +1 @@ -{ - "N_mic": 4, - "N_edge_sub": 3, - "N_edge_prod": 3, - "N_edge": 3, - "N_unbalanced": 2, - "N_enzyme": 3, - "N_er": 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_er": [ - 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": [], - "N_dgf_fixed": 1, - "dgf_fixed": [ - -1.0 - ], - "ix_dgf_free": [ - 2 - ], - "ix_dgf_fixed": [ - 1 - ], - "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, - "penalize_non_steady": 0, - "steady_state_threshold_abs": 1e-06, - "steady_state_threshold_rel": 0.001, - "steady_state_penalty_rel": 1e-08, - "rel_tol_ode": 1e-06, - "abs_tol_ode": 1e-04, - "max_num_steps_ode": 1000000, - "rel_tol_alg": 1e-06, - "abs_tol_alg": 1e-04, - "max_num_steps_alg": 1000000, - "conc_init": [ - [ - 0.59, - 0.38 - ], - [ - 0.54, - 0.38 - ] - ] -} +{"N_mic": 4, "N_pool": 0, "N_edge_sub": 3, "N_edge_prod": 3, "N_edge": 3, "N_unbalanced": 2, "N_independent": 2, "N_dependent": 0, "N_enzyme": 3, "N_er": 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]], "left_nullspace_independent": [[]], "N_reaction": 3, "N_metabolite": 2, "independent_bal_ix": [2, 3], "dependent_bal_ix": [], "unbalanced_mic_ix": [1, 4], "ci_mic_ix": [2], "edge_type": [1, 2, 1], "edge_to_enzyme": [1, 2, 3], "edge_to_er": [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": [], "N_dgf_fixed": 1, "dgf_fixed": [-1.0], "ix_dgf_free": [2], "ix_dgf_fixed": [1], "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, "penalize_non_steady": 0, "steady_state_threshold_abs": 1e-06, "steady_state_threshold_rel": 0.001, "steady_state_penalty_abs": 1e-12, "rel_tol_ode": 1e-06, "abs_tol_ode": 0.0001, "timepoint": 10000.0, "max_num_steps_ode": 1000000, "rel_tol_alg": 1e-06, "abs_tol_alg": 0.0001, "max_num_steps_alg": 1000000, "conc_init": [[0.59, 0.38], [0.54, 0.38]]} \ No newline at end of file diff --git a/maud/data_model/experiment.py b/maud/data_model/experiment.py index 843ba641c..b5f171361 100644 --- a/maud/data_model/experiment.py +++ b/maud/data_model/experiment.py @@ -12,6 +12,7 @@ class MeasurementType(str, Enum): """Possible types of measurement.""" MIC = "mic" + CONSERVED_MOIETY = "conserved_moiety" FLUX = "flux" ENZYME = "enzyme" @@ -27,6 +28,8 @@ class Measurement(BaseModel): compartment: Optional[str] = None reaction: Optional[str] = None enzyme: Optional[str] = None + moiety_id: Optional[str] = None + moiety_group: Optional[List[str]] = None @computed_field def target_id(self) -> str: @@ -38,9 +41,13 @@ def target_id(self) -> str: elif self.target_type == MeasurementType.FLUX: assert self.reaction is not None return self.reaction - else: + elif self.target_type == MeasurementType.ENZYME: assert self.enzyme is not None return self.enzyme + else: + assert self.moiety_group is not None + assert self.moiety_id is not None + return self.moiety_id class EnzymeKnockout(BaseModel): diff --git a/maud/data_model/id_component.py b/maud/data_model/id_component.py index a80ee4872..0a99b0f5d 100644 --- a/maud/data_model/id_component.py +++ b/maud/data_model/id_component.py @@ -13,3 +13,4 @@ class IdComponent(str, Enum): EXPERIMENT = "experiment" PHOSPHORYLATION_MODIFYING_ENZYME = "phosphorylation_modifying_enzyme" MODIFICATION_TYPE = "modification_type" + CONSERVED_MOIETY = "conserved_moiety" diff --git a/maud/data_model/kinetic_model.py b/maud/data_model/kinetic_model.py index f7de22f60..9e8f02e6c 100644 --- a/maud/data_model/kinetic_model.py +++ b/maud/data_model/kinetic_model.py @@ -3,6 +3,7 @@ from enum import Enum from typing import Dict, List, Optional, Union +import numpy as np import pandas as pd from pydantic import ( BaseModel, @@ -13,6 +14,7 @@ ) from maud.data_model.hardcoding import ID_SEPARATOR +from maud.utility_functions import get_left_nullspace class ReactionMechanism(int, Enum): @@ -80,6 +82,25 @@ def id_must_not_contain_seps(cls, v): return v +class ConservedMoiety(BaseModel): + """Maud representation of conserved moiety groups.""" + + id: str + name: Optional[str] + moiety_group: List[str] + + @field_validator("id") + def id_must_not_contain_seps(cls, v): + """Check that the id doesn't contain ID_SEPARATOR.""" + assert ID_SEPARATOR not in v + return v + + @computed_field + def solve(self) -> str: + """Return dependent moiety.""" + return self.moiety_group[-1] + + class Compartment(BaseModel): """Maud representation of an intra-cellular compartment. @@ -248,6 +269,7 @@ class KineticModel(BaseModel): allosteries: Optional[List[Allostery]] allosteric_enzymes: Optional[List[Enzyme]] competitive_inhibitions: Optional[List[CompetitiveInhibition]] + conserved_moiety: Optional[List[ConservedMoiety]] phosphorylations: Optional[List[Phosphorylation]] model_config: ConfigDict = {"arbitrary_types_allowed": True} @@ -268,6 +290,34 @@ def stoichiometric_matrix(self) -> pd.DataFrame: """Add the stoichiometric_matrix field.""" return get_stoichiometric_matrix(self.edges, self.mics, self.reactions) + @computed_field + def left_nullspace(self) -> pd.DataFrame: + """Calculate the left nullspace.""" + balanced_mics = [mic.id for mic in self.mics if mic.balanced] + left_nullspace = get_left_nullspace( + self.stoichiometric_matrix.loc[balanced_mics] + ) + left_nullspace[np.abs(left_nullspace) < 1e-10] = 0 + left_nullspace[np.abs(left_nullspace) > 1e-10] = 1 + + return pd.DataFrame(left_nullspace.T, columns=balanced_mics) + + @computed_field + def left_nullspace_independent(self) -> pd.DataFrame: + """Add the independent left_nullspace entries.""" + lns = self.left_nullspace + if len(lns) > 0: + dependent_mets = [mic.solve for mic in self.conserved_moiety] + dep_idx = lns[dependent_mets] + dep_idx = [ + dep_idx.index[dep_idx[met] == 1].item() + for met in dependent_mets + ] + dependent_mets = [dependent_mets[i] for i in dep_idx] + return self.left_nullspace.drop(dependent_mets, axis=1) + else: + return None + @computed_field def phosphorylation_modifying_enzymes( self, @@ -387,6 +437,32 @@ def phosphorylation_references_must_exist(self): ), f"{p.id} has bad enzyme_id" return self + @model_validator(mode="after") + def correct_conserved_moieties_defined(self): + """Compare the defined conserved_moieties to the left nullspace.""" + balanced_metabolites = [m.id for m in self.mics if m.balanced] + left_nullspace = self.left_nullspace.values + left_nullspace = left_nullspace.astype(int) + if len(left_nullspace) > 0: + mgs = [ + sorted(con_moi.moiety_group) + for con_moi in self.conserved_moiety + ] + cms = [ + sorted([m for m, n in zip(balanced_metabolites, lns) if n]) + for lns in left_nullspace + ] + for cm in cms: + assert ( + cm in mgs + ), f"{cm} not defined in kinetic model `moiety_group`" + for mg in mgs: + assert ( + mg in cms + ), f"{mg} defined in kinetic model `moiety_group` \ + but may not be a conserved moiety" + return self + def get_stoichiometric_matrix( edges: List[Union[Reaction, EnzymeReaction]], diff --git a/maud/data_model/maud_config.py b/maud/data_model/maud_config.py index b578493d5..9cd257126 100644 --- a/maud/data_model/maud_config.py +++ b/maud/data_model/maud_config.py @@ -11,6 +11,7 @@ class ODESolverConfig(BaseModel): rel_tol: float = 1e-9 abs_tol: float = 1e-9 max_num_steps: int = int(1e7) + timepoint: float = 1e4 model_config: ConfigDict = {"frozen": True} @@ -44,7 +45,7 @@ class MaudConfig(BaseModel): :param user_inits_file: path to a csv file of initial values. :param steady_state_threshold_abs: abs threshold for Sv=0 be at steady state :param steady_state_threshold_rel: rel threshold for Sv=0 be at steady state - :param steady_state_penalty_rel: relative standard deviation for SV checks when penalize_non_steady is true. + :param steady_state_penalty_abs: relative standard deviation for SV checks when penalize_non_steady is true. :param default_initial_concentration: in molecule_unit per volume_unit :param drain_small_conc_corrector: number for correcting small conc drains :param molecule_unit: A unit for counting molecules, like 'mol' or 'mmol' @@ -73,7 +74,7 @@ class MaudConfig(BaseModel): penalize_non_steady: bool = False steady_state_threshold_abs: float = 1e-8 steady_state_threshold_rel: float = 1e-3 - steady_state_penalty_rel: float = 1e-8 + steady_state_penalty_abs: float = 1e-12 default_initial_concentration: float = 0.01 drain_small_conc_corrector: float = 1e-6 molecule_unit: str = "mmol" diff --git a/maud/data_model/maud_init.py b/maud/data_model/maud_init.py index c7e2556d2..d7d6603f6 100644 --- a/maud/data_model/maud_init.py +++ b/maud/data_model/maud_init.py @@ -42,6 +42,7 @@ class InitInput(BaseModel): drain: ParamInitInput = None conc_enzyme: ParamInitInput = None conc_pme: ParamInitInput = None + conc_moiety_pool: ParamInitInput = None def get_init_atom_input_ids( diff --git a/maud/data_model/maud_parameter.py b/maud/data_model/maud_parameter.py index 2f2ad0fb5..86911bd6a 100644 --- a/maud/data_model/maud_parameter.py +++ b/maud/data_model/maud_parameter.py @@ -353,6 +353,20 @@ class ConcUnbalanced(TrainTestParameter): default_scale: float = 2.0 +class ConservedMoiety(TrainTestParameter): + """Parent class for unbalanced mic concentration parameters.""" + + name: str = "conc_moiety_pool_train" + shape_names: List[str] = ["N_experiment_train", "N_moiety_pool"] + id_components: List[List[IdComponent]] = [ + [IdComponent.EXPERIMENT], + [IdComponent.CONSERVED_MOIETY], + ] + non_negative: bool = True + default_loc: float = -2.3 + default_scale: float = 2.0 + + class ConcPme(TrainTestParameter): """Parent class for pme concentration parameters.""" diff --git a/maud/data_model/parameter_input.py b/maud/data_model/parameter_input.py index 7e144c04b..ae2711935 100644 --- a/maud/data_model/parameter_input.py +++ b/maud/data_model/parameter_input.py @@ -21,6 +21,7 @@ class ParameterInputAtom(BaseModel): experiment: Optional[str] = None phosphorylation_modifying_enzyme: Optional[str] = None modification_type: Optional[str] = None + conserved_moiety: Optional[str] = None location: Optional[float] = None exploc: Optional[float] = None scale: Optional[PositiveFloat] = None @@ -96,3 +97,6 @@ class ParameterSetInput(BaseModel): default_factory=list ) conc_pme: Optional[List[ParameterInputAtom]] = Field(default_factory=list) + conc_moiety_pool: Optional[List[ParameterInputAtom]] = Field( + default_factory=list + ) diff --git a/maud/data_model/parameter_set.py b/maud/data_model/parameter_set.py index cc207c75c..21cdb6578 100644 --- a/maud/data_model/parameter_set.py +++ b/maud/data_model/parameter_set.py @@ -226,6 +226,32 @@ def conc_enzyme_test(self) -> mp.ConcEnzyme: """Add the conc_enzyme_test field.""" return self._get_conc_enzyme(train=False) + def _get_conc_moiety_pool(self, train: bool) -> mp.ConservedMoiety: + exp_ids = self._get_experiments(train) + if len(self.kinetic_model.left_nullspace) > 0: + conserved_moiety_ids = [ + cm.id for cm in self.kinetic_model.conserved_moiety + ] + result = mp.ConservedMoiety( + ids=[exp_ids, conserved_moiety_ids], + split_ids=[[exp_ids], [conserved_moiety_ids]], + user_input=self.parameter_set_input.conc_moiety_pool, + init_input=self.init_input.conc_moiety_pool, + ) + return result if train else result.test() + else: + return [] + + @computed_field + def conc_moiety_pool_train(self) -> mp.ConservedMoiety: + """Add the conc_enzyme_train field.""" + return self._get_conc_moiety_pool(train=True) + + @computed_field + def conc_moiety_pool_test(self) -> mp.ConservedMoiety: + """Add the conc_enzyme_test field.""" + return self._get_conc_moiety_pool(train=False) + def _get_conc_unbalanced(self, train: bool) -> mp.ConcUnbalanced: exp_ids = self._get_experiments(train) measurements = self._get_measurements(train, MeasurementType.MIC) diff --git a/maud/getting_idatas.py b/maud/getting_idatas.py index dbb5aa674..0c8b7a486 100644 --- a/maud/getting_idatas.py +++ b/maud/getting_idatas.py @@ -29,6 +29,11 @@ def get_idata(csvs: List[str], mi: MaudInput, mode: str) -> az.InferenceData: MeasurementType.ENZYME, ] ) + dependent_mics = ( + [m.solve for m in mi.kinetic_model.conserved_moiety] + if mi.kinetic_model.conserved_moiety is not None + else [] + ) coords = { "enzymes": [e.id for e in mi.kinetic_model.enzymes], "experiments": [e.id for e in experiments], @@ -36,12 +41,21 @@ def get_idata(csvs: List[str], mi: MaudInput, mode: str) -> az.InferenceData: "drains": [r.id for r in mi.kinetic_model.drains], "metabolites": [m.id for m in mi.kinetic_model.metabolites], "mics": [m.id for m in mi.kinetic_model.mics], + "conserved_moieties": ( + [cm.id for cm in mi.kinetic_model.conserved_moiety] + if mi.kinetic_model.conserved_moiety is not None + else [] + ), "edges": [e.id for e in mi.kinetic_model.edges], "edges1": [e.id for e in mi.kinetic_model.edges], "unbalanced_mics": [ m.id for m in mi.kinetic_model.mics if not m.balanced ], - "balanced_mics": [m.id for m in mi.kinetic_model.mics if m.balanced], + "independent_mics": [ + m.id + for m in mi.kinetic_model.mics + if (m.balanced) & (m.id not in dependent_mics) + ], "phosphorylations": ( [p.id for p in mi.kinetic_model.phosphorylations] if mi.kinetic_model.phosphorylations is not None @@ -80,6 +94,7 @@ def get_idata(csvs: List[str], mi: MaudInput, mode: str) -> az.InferenceData: dims = { f"flux_{mode}": ["experiments", "reactions"], f"conc_{mode}": ["experiments", "mics"], + f"conc_moiety_pool_{mode}": ["experiments", "conserved_moieties"], f"log_conc_enzyme_{mode}_z": ["experiments", "enzymes"], f"conc_enzyme_{mode}": ["experiments", "enzymes"], f"conc_unbalanced_{mode}": ["experiments", "unbalanced_mics"], @@ -109,14 +124,14 @@ def get_idata(csvs: List[str], mi: MaudInput, mode: str) -> az.InferenceData: f"llik_flux_{mode}": ["yfluxs"], "concentration_control_matrix": [ "experiments", - "balanced_mics", + "independent_mics", "edges", ], "flux_control_matrix": ["experiments", "edges", "edges1"], "flux_response_coefficient": ["experiments", "edges", "enzymes"], "concentration_response_coefficient": [ "experiments", - "balanced_mics", + "independent_mics", "enzymes", ], } diff --git a/maud/getting_stan_inputs.py b/maud/getting_stan_inputs.py index 3763c1e21..4469d6832 100644 --- a/maud/getting_stan_inputs.py +++ b/maud/getting_stan_inputs.py @@ -2,6 +2,7 @@ from typing import Dict, Iterable, List, Tuple, Union +import pandas as pd from scipy.stats import gmean from maud.data_model.experiment import Experiment, MeasurementType @@ -9,6 +10,7 @@ from maud.data_model.kinetic_model import ( Allostery, CompetitiveInhibition, + ConservedMoiety, Enzyme, EnzymeReaction, KineticModel, @@ -136,6 +138,18 @@ def get_network_properties_input( er_codes = codify_maud_object(kinetic_model.ers) km_ids = parameters.km.ids[0] km_codes = dict(zip(km_ids, range(1, len(km_ids) + 1))) + if len(kinetic_model.left_nullspace) > 0: + lns_ind = kinetic_model.left_nullspace_independent.values.tolist() + dependent_mets = [mic.solve for mic in kinetic_model.conserved_moiety] + dep_idx = kinetic_model.left_nullspace[dependent_mets] + dep_idx = [ + dep_idx.index[dep_idx[met] == 1].item() for met in dependent_mets + ] + dependent_mets = [dependent_mets[i] for i in dep_idx] + dependent_mic_codes = [mic_codes[lns] for lns in dependent_mets] + else: + lns_ind = [[]] + dependent_mic_codes = [] if allosteries is not None: tc_codes = codify_maud_object( [ @@ -209,7 +223,11 @@ def get_network_properties_input( for e in edges ] mic_met_code = [metabolite_codes[m.metabolite_id] for m in mics] - balanced_mic_codes = [mic_codes[mic.id] for mic in mics if mic.balanced] + independent_mic_codes = [ + mic_codes[mic.id] + for mic in mics + if ((mic.balanced) & (mic_codes[mic.id] not in dependent_mic_codes)) + ] unbalanced_mic_codes = [ mic_codes[mic.id] for mic in mics if not mic.balanced ] @@ -337,10 +355,13 @@ def get_network_properties_input( ) return { "N_mic": len(kinetic_model.mics), + "N_pool": len(dependent_mic_codes), "N_edge_sub": len(sub_by_edge_long), "N_edge_prod": len(prod_by_edge_long), "N_edge": len(S.columns), "N_unbalanced": len(unbalanced_mic_codes), + "N_independent": len(independent_mic_codes), + "N_dependent": len(dependent_mic_codes), "N_enzyme": len(enzymes), "N_er": len(kinetic_model.ers), "N_phosphorylation": len(phosphorylation_codes), @@ -353,9 +374,11 @@ def get_network_properties_input( "N_sub_km": len(sub_km_ix_by_edge_long), "N_prod_km": len(prod_km_ix_by_edge_long), "S": S.values.tolist(), + "left_nullspace_independent": lns_ind, "N_reaction": len(reactions), "N_metabolite": len(metabolite_codes), - "balanced_mic_ix": balanced_mic_codes, + "independent_bal_ix": independent_mic_codes, + "dependent_bal_ix": dependent_mic_codes, "unbalanced_mic_ix": unbalanced_mic_codes, "ci_mic_ix": ci_mic_codes, "edge_type": edge_mechanism_code, @@ -518,9 +541,10 @@ def get_config_input(config: MaudConfig): "penalize_non_steady": int(config.penalize_non_steady), "steady_state_threshold_abs": config.steady_state_threshold_abs, "steady_state_threshold_rel": config.steady_state_threshold_rel, - "steady_state_penalty_rel": config.steady_state_penalty_rel, + "steady_state_penalty_abs": config.steady_state_penalty_abs, "rel_tol_ode": config.ode_solver_config.rel_tol, "abs_tol_ode": config.ode_solver_config.abs_tol, + "timepoint": config.ode_solver_config.timepoint, "max_num_steps_ode": config.ode_solver_config.max_num_steps, "rel_tol_alg": config.algebra_solver_config.rel_tol, "abs_tol_alg": config.algebra_solver_config.abs_tol, @@ -559,22 +583,41 @@ def get_conc_init( """ default = config.default_initial_concentration balanced_mics = [m for m in kinetic_model.mics if m.balanced] - inits = [[default for _ in balanced_mics] for e in experiments] - for i, experiment in enumerate(experiments): + balanced_mic_ids = [m.id for m in balanced_mics] + experiment_ids = [exp.id for exp in experiments] + experiments_train = [exp.id for exp in experiments if exp.is_train] + experiments_test = [exp.id for exp in experiments if exp.is_test] + if len(kinetic_model.left_nullspace) > 0: + dependent_mics = [m.solve for m in kinetic_model.conserved_moiety] + moiety_groups = [ + mg.moiety_group for mg in kinetic_model.conserved_moiety + ] + else: + dependent_mics = [] + moiety_groups = [] + independent_mics = [ + m.id + for m in kinetic_model.mics + if (m.balanced) & (m.id not in dependent_mics) + ] + inits = pd.DataFrame( + default, index=experiment_ids, columns=balanced_mic_ids + ) + for experiment in experiments: msmts = [ m for m in experiment.measurements if m.target_type == MeasurementType.MIC ] + experiment.initial_state - for j, mic in enumerate(balanced_mics): + for mic in balanced_mics: if any(mic.id == m.target_id for m in msmts): - inits[i][j] = gmean( + inits.loc[experiment.id, mic.id] = gmean( [m.value for m in msmts if mic.id == m.target_id] ) - inits_train = [ - inits[i] for i, exp in enumerate(experiments) if exp.is_train - ] - inits_test = [inits[i] for i, exp in enumerate(experiments) if exp.is_test] + for mg in moiety_groups: + inits[mg] = inits[mg].apply(lambda x: x / x.sum(), axis=1) + inits_train = inits.loc[experiments_train, independent_mics].values + inits_test = inits.loc[experiments_test, independent_mics].values return inits_train, inits_test @@ -589,6 +632,7 @@ def get_conc_init( Reaction, Experiment, CompetitiveInhibition, + ConservedMoiety, ] diff --git a/maud/parsing_kinetic_models.py b/maud/parsing_kinetic_models.py index 9902c40e7..e9ab46fcd 100644 --- a/maud/parsing_kinetic_models.py +++ b/maud/parsing_kinetic_models.py @@ -7,6 +7,7 @@ Allostery, Compartment, CompetitiveInhibition, + ConservedMoiety, Enzyme, EnzymeReaction, KineticModel, @@ -87,6 +88,17 @@ def parse_kinetic_model(raw: dict) -> KineticModel: ] else: competitive_inhibitions = None + if "conserved_moiety" in raw.keys(): + conserved_moiety = [ + ConservedMoiety( + id=cm["id"], + name=read_with_fallback("name", cm, None), + moiety_group=cm["moiety_group"], + ) + for cm in raw["conserved_moiety"] + ] + else: + conserved_moiety = None return KineticModel( name=name, metabolites=metabolites, @@ -98,6 +110,7 @@ def parse_kinetic_model(raw: dict) -> KineticModel: allosteries=allosteries, allosteric_enzymes=allosteric_enzymes, competitive_inhibitions=competitive_inhibitions, + conserved_moiety=conserved_moiety, phosphorylations=phosphorylations, ) diff --git a/maud/running_stan.py b/maud/running_stan.py index e61285101..cfdef7599 100644 --- a/maud/running_stan.py +++ b/maud/running_stan.py @@ -311,13 +311,13 @@ def predict( "output_dir": output_dir, "iter_warmup": 0, "iter_sampling": 1, + "adapt_engaged": False, "fixed_param": True, - "show_progress": False, } if mi.config.cmdstanpy_config_predict is not None: sample_args = { - **sample_args, **mi.config.cmdstanpy_config_predict, + **sample_args, } sample_args["chains"] = 1 mcmc_draw = model.sample(**sample_args) diff --git a/maud/stan/functions.stan b/maud/stan/functions.stan index cf7921834..96e93e283 100644 --- a/maud/stan/functions.stan +++ b/maud/stan/functions.stan @@ -333,8 +333,7 @@ functions { array[] int phosphorylation_pme) { int N_edge = cols(S); vector[N_edge] vmax = get_vmax_by_edge(enzyme, kcat, edge_to_enzyme, - edge_to_er, - edge_type); + edge_to_er, edge_type); vector[N_edge] reversibility = get_reversibility(dgr, temperature, S, conc, edge_type); vector[N_edge] free_enzyme_ratio = get_free_enzyme_ratio(conc, S, km, ki, @@ -376,12 +375,13 @@ functions { .* phosphorylation .* drain_by_edge; } - vector dbalanced_dt(real time, vector current_balanced, vector unbalanced, - array[] int balanced_ix, array[] int unbalanced_ix, + vector dbalanced_dt(real time, vector current_independent, vector unbalanced, + vector conc_pool, array[] int independent_bal_ix, array[] int dependent_bal_ix, + array[] int unbalanced_ix, vector enzyme, vector dgr, vector kcat, vector km, vector ki, vector tc, vector dc, vector kcat_pme, vector conc_pme, vector drain, real temperature, - real drain_small_conc_corrector, matrix S, + real drain_small_conc_corrector, matrix S, matrix left_nullspace_independent, vector subunits, array[] int edge_type, array[] int edge_to_enzyme, array[] int edge_to_er, array[] int edge_to_drain, array[] int ci_mic_ix, @@ -402,8 +402,9 @@ functions { array[,] int phosphorylation_ix_bounds, array[] int phosphorylation_type, array[] int phosphorylation_pme) { - vector[rows(current_balanced) + rows(unbalanced)] current_concentration; - current_concentration[balanced_ix] = current_balanced; + vector[size(independent_bal_ix) + size(dependent_bal_ix)+ rows(unbalanced)] current_concentration; + current_concentration[independent_bal_ix] = current_independent; + current_concentration[dependent_bal_ix] = conc_pool - left_nullspace_independent * current_independent; current_concentration[unbalanced_ix] = unbalanced; vector[cols(S)] edge_flux = get_edge_flux(current_concentration, enzyme, dgr, kcat, km, ki, tc, dc, @@ -430,7 +431,7 @@ functions { phosphorylation_ix_bounds, phosphorylation_type, phosphorylation_pme); - return (S * edge_flux)[balanced_ix]; + return S[independent_bal_ix] * edge_flux; } complex_vector get_complex_edge_flux_enzyme(vector conc, complex_vector enzyme, @@ -530,19 +531,19 @@ functions { return get_imag(edge_flux) ./ get_imag(complex_enzyme_step); } matrix get_concentration_control_matrix(matrix S, matrix elasticity, - array[] int balanced_mic_ix, - int N_edge, int N_balanced) { - matrix[N_balanced, N_edge] balanced_S = S[balanced_mic_ix, : ]; - return -generalized_inverse(balanced_S * elasticity) * balanced_S; + array[] int independent_bal_ix, + int N_edge, int N_independent) { + matrix[N_independent, N_edge] independent_bal_S= S[independent_bal_ix, : ]; + return -generalized_inverse(independent_bal_S * elasticity) * independent_bal_S; } matrix get_flux_control_matrix(matrix S, matrix elasticity, - array[] int balanced_mic_ix, int N_edge, - int N_balanced) { - matrix[N_balanced, N_edge] balanced_S = S[balanced_mic_ix, : ]; + array[] int independent_bal_ix, int N_edge, + int N_independent) { + matrix[N_independent, N_edge] independent_bal_S= S[independent_bal_ix, : ]; matrix[N_edge, N_edge] I = identity_matrix(N_edge); return I - - elasticity * generalized_inverse(balanced_S * elasticity) - * balanced_S; + - elasticity * generalized_inverse(independent_bal_S * elasticity) + * independent_bal_S; } complex_vector get_complex_edge_flux_metabolite(complex_vector conc, vector enzyme, vector dgr, @@ -579,7 +580,7 @@ functions { array[] int phosphorylation_type, array[] int phosphorylation_pme) { int N_edge = cols(S); - vector[N_edge] vmax = get_vmax_by_edge(enzyme, kcat, edge_to_enzyme,edge_to_er, + vector[N_edge] vmax = get_vmax_by_edge(enzyme, kcat, edge_to_enzyme, edge_to_er, edge_type); complex_vector[N_edge] reversibility = get_reversibility_complex(dgr, temperature, @@ -805,10 +806,12 @@ functions { } return out; } - -vector maud_ae_system(vector conc, +/* +vector maud_ae_system(vector conc_ind, data real rel_tol_ode, data real abs_tol_ode, data int max_num_steps_ode, - vector conc_unbalanced, + vector init_dc_ind_dt, + vector conc_unbalanced, vector pool, array[] int independent_bal_ix, + array[] int dependent_bal_ix, array[] int balanced_mic_ix, array[] int unbalanced_mic_ix, vector conc_enzyme, vector dgr, vector kcat, vector km, vector ki, vector tc, vector dc, vector kcat_pme, @@ -834,9 +837,20 @@ vector maud_ae_system(vector conc, array[,] int phosphorylation_ix_bounds, array[] int phosphorylation_type, array[] int phosphorylation_pme){ - return conc - to_vector(ode_bdf_tol(dbalanced_dt, conc, 0, {1}, + \* + Delaying implementation of AlgebraSolver + DAESolver until we can confirm + DAE solver works independently. Concern with AlgebraSolver is that the independent + concentrations may be chosen higher than the selected pool size. Potential solutions + to overcome there are to sample the concentraitons as part of a simplex. Additionally + we could always select the smallest concentration as independent so the solver is more + unlikely to sample over the bound. + *\ + return conc - to_vector(dae_tol(dae_residual, conc, init_dc_ind_dt, 0, {1}, rel_tol_ode, abs_tol_ode, max_num_steps_ode, conc_unbalanced, + pool, + independent_bal_ix, + dependent_bal_ix balanced_mic_ix, unbalanced_mic_ix, conc_enzyme, @@ -866,5 +880,6 @@ vector maud_ae_system(vector conc, phosphorylation_type, phosphorylation_pme)[1]); } +*/ } diff --git a/maud/stan/model.stan b/maud/stan/model.stan index e9af4a459..090e8605a 100644 --- a/maud/stan/model.stan +++ b/maud/stan/model.stan @@ -2,9 +2,12 @@ data { // network properties int N_mic; + int N_pool; int N_edge_sub; int N_edge_prod; int N_unbalanced; + int N_independent; + int N_dependent; int N_metabolite; int N_km; int N_sub_km; @@ -21,7 +24,9 @@ data { int N_competitive_inhibition; int N_dgf_fixed; matrix[N_mic, N_edge] S; - array[N_mic - N_unbalanced] int balanced_mic_ix; + matrix[N_pool, N_independent] left_nullspace_independent; + array[N_independent] int independent_bal_ix; + array[N_dependent] int dependent_bal_ix; array[N_unbalanced] int unbalanced_mic_ix; array[N_competitive_inhibition] int ci_mic_ix; array[N_edge] int edge_type; // 1 = reversible modular rate law, 2 = drain @@ -94,8 +99,10 @@ data { array[2, N_experiment_train] vector[N_unbalanced] priors_conc_unbalanced_train; array[2, N_experiment_train] vector[N_enzyme] priors_conc_enzyme_train; array[2, N_experiment_train] vector[N_drain] priors_drain_train; + array[2, N_experiment_train] vector[N_pool] priors_conc_moiety_pool_train; // configuration - array[N_experiment_train] vector[N_mic - N_unbalanced] conc_init; + real timepoint; + array[N_experiment_train] vector[N_independent] conc_init; real rel_tol_ode; real abs_tol_ode; int max_num_steps_ode; @@ -104,7 +111,7 @@ data { int max_num_steps_alg; real steady_state_threshold_abs; real steady_state_threshold_rel; - real steady_state_penalty_rel; + real steady_state_penalty_abs; int likelihood; // set to 0 for priors-only mode real drain_small_conc_corrector; int penalize_non_steady; @@ -119,7 +126,6 @@ transformed data { } matrix[N_metabolite-N_dgf_fixed, N_metabolite-N_dgf_fixed] prior_cov_dgf_free_chol = cholesky_decompose(prior_cov_dgf_free); - int N_balanced = N_mic - N_unbalanced; } parameters { vector[N_metabolite-N_dgf_fixed] dgf_free; @@ -134,6 +140,7 @@ parameters { array[N_experiment_train] vector[N_enzyme] log_conc_enzyme_train_z; array[N_experiment_train] vector[N_pme] log_conc_pme_train_z; array[N_experiment_train] vector[N_unbalanced] log_conc_unbalanced_train_z; + array[N_experiment_train] vector[N_pool] log_conc_moiety_pool_train_z; } transformed parameters { // unfix @@ -157,19 +164,25 @@ transformed parameters { log_conc_unbalanced_train_z); array[N_experiment_train] vector[N_pme] conc_pme_train = unz_log_2d(priors_conc_pme_train, log_conc_pme_train_z); + array[N_experiment_train] vector[N_pool] conc_moiety_pool_train = unz_log_2d(priors_conc_moiety_pool_train, + log_conc_moiety_pool_train_z); // transform array[N_experiment_train] vector[N_mic] conc_train; array[N_experiment_train] vector[N_reaction] flux_train; array[N_experiment_train] vector[N_edge] dgr_train; - matrix[N_experiment_train, N_mic - N_unbalanced] steady_dev; + matrix[N_experiment_train, N_independent] steady_dev; for (e in 1 : N_experiment_train) { + vector[N_independent] conc_init_exp = conc_init[e]; dgr_train[e] = get_dgr(S, dgf, temperature_train[e], mic_to_met, water_stoichiometry, transported_charge, psi_train[e]); flux_train[e] = rep_vector(0, N_reaction); vector[N_enzyme] conc_enzyme_experiment = conc_enzyme_train[e]; + vector[N_unbalanced] conc_unbalanced_experiment = conc_unbalanced_train[e]; + vector[N_drain] drain_experiment = drain_train[e]; + vector[N_pool] conc_moiety_pool_experiment = conc_moiety_pool_train[e]; vector[N_pme] conc_pme_experiment = conc_pme_train[e]; - vector[N_mic - N_unbalanced] conc_balanced_experiment; + vector[N_independent] conc_independent_balanced_experiment; int N_eko_experiment = measure_ragged(enzyme_knockout_train_bounds, e); int N_pko_experiment = measure_ragged(pme_knockout_train_bounds, e); if (N_eko_experiment > 0) { @@ -185,43 +198,51 @@ transformed parameters { e); conc_pme_experiment[pko_experiment] = rep_vector(0, N_pko_experiment); } - conc_balanced_experiment = solve_newton_tol(maud_ae_system, - conc_init[e], - rel_tol_alg, abs_tol_alg, - max_num_steps_alg, - rel_tol_ode, abs_tol_ode, - max_num_steps_ode, - conc_unbalanced_train[e], - balanced_mic_ix, - unbalanced_mic_ix, - conc_enzyme_experiment, - dgr_train[e], kcat, km, ki, - transfer_constant, - dissociation_constant, kcat_pme, - conc_pme_experiment, - drain_train[e], - temperature_train[e], - drain_small_conc_corrector, S, - subunits, edge_type, - edge_to_enzyme, edge_to_er, edge_to_drain, - ci_mic_ix, sub_km_ix_by_edge_long, - sub_km_ix_by_edge_bounds, - prod_km_ix_by_edge_long, - prod_km_ix_by_edge_bounds, - sub_by_edge_long, - sub_by_edge_bounds, - prod_by_edge_long, - prod_by_edge_bounds, ci_ix_long, - ci_ix_bounds, allostery_ix_long, - allostery_ix_bounds, - allostery_type, allostery_mic, - edge_to_tc, - phosphorylation_ix_long, - phosphorylation_ix_bounds, - phosphorylation_type, - phosphorylation_pme); - conc_train[e, balanced_mic_ix] = conc_balanced_experiment; - conc_train[e, unbalanced_mic_ix] = conc_unbalanced_train[e]; + for (p in 1 : N_pool){ + for (m in 1 : N_independent){ + if (left_nullspace_independent[p,m] == 1){ + conc_init_exp[m] = conc_init_exp[m] * conc_moiety_pool_experiment[p]; + } + } + } + conc_independent_balanced_experiment = ode_bdf_tol(dbalanced_dt, conc_init_exp, 0, {timepoint}, + rel_tol_ode, abs_tol_ode, max_num_steps_ode, + conc_unbalanced_experiment, + conc_moiety_pool_experiment, + independent_bal_ix, + dependent_bal_ix, + unbalanced_mic_ix, + conc_enzyme_experiment, + dgr_train[e], kcat, km, ki, + transfer_constant, + dissociation_constant, kcat_pme, + conc_pme_experiment, + drain_experiment, + temperature_train[e], + drain_small_conc_corrector, S, + left_nullspace_independent, + subunits, edge_type, + edge_to_enzyme, edge_to_er, edge_to_drain, + ci_mic_ix, sub_km_ix_by_edge_long, + sub_km_ix_by_edge_bounds, + prod_km_ix_by_edge_long, + prod_km_ix_by_edge_bounds, + sub_by_edge_long, + sub_by_edge_bounds, + prod_by_edge_long, + prod_by_edge_bounds, ci_ix_long, + ci_ix_bounds, allostery_ix_long, + allostery_ix_bounds, + allostery_type, allostery_mic, + edge_to_tc, + phosphorylation_ix_long, + phosphorylation_ix_bounds, + phosphorylation_type, + phosphorylation_pme)[1]; + conc_train[e, independent_bal_ix] = conc_independent_balanced_experiment; + conc_train[e, dependent_bal_ix] = conc_moiety_pool_experiment - left_nullspace_independent * conc_independent_balanced_experiment; + conc_train[e, unbalanced_mic_ix] = conc_unbalanced_experiment; + { vector[N_edge] edge_flux = get_edge_flux(conc_train[e], conc_enzyme_experiment, @@ -252,7 +273,7 @@ transformed parameters { phosphorylation_ix_bounds, phosphorylation_type, phosphorylation_pme); - steady_dev[e] = (S * edge_flux)[balanced_mic_ix]'; + steady_dev[e] = (S * edge_flux)[independent_bal_ix]'; for (j in 1 : N_edge) { flux_train[e, edge_to_reaction[j]] += edge_flux[j]; } @@ -272,6 +293,7 @@ model { log_conc_unbalanced_train_z[ex] ~ std_normal(); log_conc_enzyme_train_z[ex] ~ std_normal(); log_conc_pme_train_z[ex] ~ std_normal(); + log_conc_moiety_pool_train_z[ex] ~ std_normal(); drain_train_z[ex] ~ std_normal(); psi_train_z[ex] ~ std_normal(); } @@ -290,9 +312,7 @@ model { } if (penalize_non_steady == 1) { for (xpt in 1 : N_experiment_train) { - steady_dev[xpt] ~ normal(0.0, - conc_train[xpt, balanced_mic_ix] - * steady_state_penalty_rel); + 0 ~ normal(steady_dev[xpt], steady_state_penalty_abs); } } } @@ -302,12 +322,12 @@ generated quantities { vector[N_flux_measurement_train] yrep_flux_train; vector[N_conc_measurement_train] llik_conc_train; vector[N_flux_measurement_train] llik_flux_train; - array[N_experiment_train] matrix[N_mic - N_unbalanced, N_edge] concentration_control_matrix; + array[N_experiment_train] matrix[N_independent, N_edge] concentration_control_matrix; array[N_experiment_train] matrix[N_edge, N_edge] flux_control_matrix; - array[N_experiment_train] matrix[N_edge, N_balanced] elasticity; + array[N_experiment_train] matrix[N_edge, N_independent] elasticity; array[N_experiment_train] matrix[N_edge, N_enzyme] sensitivity; array[N_experiment_train] matrix[N_edge, N_enzyme] flux_response_coefficient; - array[N_experiment_train] matrix[N_balanced, N_enzyme] concentration_response_coefficient; + array[N_experiment_train] matrix[N_independent, N_enzyme] concentration_response_coefficient; array[N_experiment_train] vector[N_edge] free_enzyme_ratio_train; array[N_experiment_train] vector[N_edge] saturation_train; array[N_experiment_train] vector[N_edge] allostery_train; @@ -370,8 +390,8 @@ generated quantities { // Calculating control coefficients for (e in 1 : N_experiment_train) { vector[N_pme] conc_pme_experiment = conc_pme_train[e]; - for (bal_met in 1 : N_balanced) { - int met_idx = balanced_mic_ix[bal_met]; + for (ind_met in 1 : N_independent) { + int met_idx = independent_bal_ix[ind_met]; complex_vector[N_mic] complex_step_conc = conc_train[e]; complex_step_conc[met_idx] = complex_step_conc[met_idx] + complex_step; complex_vector[N_edge] edge_flux_met = get_complex_edge_flux_metabolite(complex_step_conc, @@ -412,7 +432,7 @@ generated quantities { phosphorylation_ix_bounds, phosphorylation_type, phosphorylation_pme); - elasticity[e][ : , bal_met] = get_flux_jacobian(edge_flux_met, + elasticity[e][ : , ind_met] = get_flux_jacobian(edge_flux_met, complex_step); } for (enz in 1 : N_enzyme) { @@ -461,12 +481,12 @@ generated quantities { } concentration_control_matrix[e] = get_concentration_control_matrix(S, elasticity[e], - balanced_mic_ix, + independent_bal_ix, N_edge, - N_balanced); + N_independent); flux_control_matrix[e] = get_flux_control_matrix(S, elasticity[e], - balanced_mic_ix, N_edge, - N_balanced); + independent_bal_ix, N_edge, + N_independent); flux_response_coefficient[e] = flux_control_matrix[e] * sensitivity[e]; concentration_response_coefficient[e] = concentration_control_matrix[e] * sensitivity[e]; diff --git a/maud/stan/out_of_sample_model.stan b/maud/stan/out_of_sample_model.stan index 713c0fa88..402872e9e 100644 --- a/maud/stan/out_of_sample_model.stan +++ b/maud/stan/out_of_sample_model.stan @@ -2,9 +2,12 @@ data { // network properties int N_mic; + int N_pool; int N_edge_sub; int N_edge_prod; int N_unbalanced; + int N_independent; + int N_dependent; int N_metabolite; int N_km; int N_sub_km; @@ -20,7 +23,9 @@ data { int N_pme; // phosphorylation modifying enzyme int N_competitive_inhibition; matrix[N_mic, N_edge] S; - array[N_mic - N_unbalanced] int balanced_mic_ix; + matrix[N_pool, N_independent] left_nullspace_independent; + array[N_independent] int independent_bal_ix; + array[N_dependent] int dependent_bal_ix; array[N_unbalanced] int unbalanced_mic_ix; array[N_competitive_inhibition] int ci_mic_ix; array[N_edge] int edge_type; // 1 = reversible modular rate law, 2 = drain @@ -60,22 +65,21 @@ data { array[N_pme_knockout_test] int pme_knockout_test_long; array[N_experiment_test, 2] int pme_knockout_test_bounds; vector[N_experiment_test] temperature_test; - // hardcoded priors + // hardcoded priors (boundary conditions) array[2, N_experiment_test] vector[N_pme] priors_conc_phos_test; array[2, N_experiment_test] vector[N_unbalanced] priors_conc_unbalanced_test; + array[2, N_experiment_test] vector[N_pool] priors_conc_moiety_pool_test; array[2, N_experiment_test] vector[N_enzyme] priors_conc_enzyme_test; array[2, N_experiment_test] vector[N_drain] priors_drain_test; // configuration - array[N_experiment_test] vector[N_mic - N_unbalanced] conc_init; + real timepoint; + array[N_experiment_test] vector[N_independent] conc_init; real rel_tol_ode; real abs_tol_ode; int max_num_steps_ode; real rel_tol_alg; real abs_tol_alg; int max_num_steps_alg; - real steady_state_threshold_abs; - real steady_state_threshold_rel; - real steady_state_penalty_rel; int likelihood; // set to 0 for priors-only mode real drain_small_conc_corrector; int penalize_non_steady; @@ -97,6 +101,7 @@ generated quantities { array[N_experiment_test] vector[N_reaction] flux_test; array[N_experiment_test] vector[N_pme] conc_pme_test; array[N_experiment_test] vector[N_unbalanced] conc_unbalanced_test; + array[N_experiment_test] vector[N_pool] conc_moiety_pool_test; array[N_experiment_test] vector[N_enzyme] conc_enzyme_test; array[N_experiment_test] vector[N_drain] drain_test; array[N_experiment_test] vector[N_edge] free_enzyme_ratio_test; @@ -112,6 +117,8 @@ generated quantities { priors_conc_phos_test[2, e])); conc_unbalanced_test[e] = to_vector(lognormal_rng(priors_conc_unbalanced_test[1, e], priors_conc_unbalanced_test[2, e])); + conc_moiety_pool_test[e] = to_vector(lognormal_rng(priors_conc_moiety_pool_test[1, e], + priors_conc_moiety_pool_test[2, e])); conc_enzyme_test[e] = to_vector(lognormal_rng(priors_conc_enzyme_test[1, e], priors_conc_enzyme_test[2, e])); } @@ -120,7 +127,7 @@ generated quantities { flux_test[e] = rep_vector(0, N_reaction); vector[N_enzyme] conc_enzyme_experiment = conc_enzyme_test[e]; vector[N_pme] conc_pme_experiment = conc_pme_test[e]; - vector[N_mic - N_unbalanced] conc_balanced_experiment; + vector[N_independent] conc_independent_balanced_experiment; int N_eko_experiment = measure_ragged(enzyme_knockout_test_bounds, e); int N_pko_experiment = measure_ragged(pme_knockout_test_bounds, e); if (N_eko_experiment > 0) { @@ -136,42 +143,42 @@ generated quantities { e); conc_pme_experiment[pko_experiment] = rep_vector(0, N_pko_experiment); } - conc_balanced_experiment = solve_newton_tol(maud_ae_system, - conc_init[e], - rel_tol_alg, abs_tol_alg, - max_num_steps_alg, - rel_tol_ode, abs_tol_ode, - max_num_steps_ode, - conc_unbalanced_test[e], - balanced_mic_ix, - unbalanced_mic_ix, - conc_enzyme_experiment, - dgr_test[e], kcat, km, ki, - transfer_constant, - dissociation_constant, kcat_pme, - conc_pme_experiment, - drain_test[e], - temperature_test[e], - drain_small_conc_corrector, S, - subunits, edge_type, - edge_to_enzyme, edge_to_er, edge_to_drain, - ci_mic_ix, sub_km_ix_by_edge_long, - sub_km_ix_by_edge_bounds, - prod_km_ix_by_edge_long, - prod_km_ix_by_edge_bounds, - sub_by_edge_long, - sub_by_edge_bounds, - prod_by_edge_long, - prod_by_edge_bounds, ci_ix_long, - ci_ix_bounds, allostery_ix_long, - allostery_ix_bounds, - allostery_type, allostery_mic, - edge_to_tc, - phosphorylation_ix_long, - phosphorylation_ix_bounds, - phosphorylation_type, - phosphorylation_pme); - conc_test[e, balanced_mic_ix] = conc_balanced_experiment; + conc_independent_balanced_experiment = ode_bdf_tol(dbalanced_dt, conc_init[e], 0, {timepoint}, + rel_tol_ode, abs_tol_ode, max_num_steps_ode, + conc_unbalanced_test[e], + conc_moiety_pool_test[e], + independent_bal_ix, + dependent_bal_ix, + unbalanced_mic_ix, + conc_enzyme_experiment, + dgr_test[e], kcat, km, ki, + transfer_constant, + dissociation_constant, kcat_pme, + conc_pme_experiment, + drain_test[e], + temperature_test[e], + drain_small_conc_corrector, S, + left_nullspace_independent, + subunits, edge_type, + edge_to_enzyme, edge_to_er, edge_to_drain, + ci_mic_ix, sub_km_ix_by_edge_long, + sub_km_ix_by_edge_bounds, + prod_km_ix_by_edge_long, + prod_km_ix_by_edge_bounds, + sub_by_edge_long, + sub_by_edge_bounds, + prod_by_edge_long, + prod_by_edge_bounds, ci_ix_long, + ci_ix_bounds, allostery_ix_long, + allostery_ix_bounds, + allostery_type, allostery_mic, + edge_to_tc, + phosphorylation_ix_long, + phosphorylation_ix_bounds, + phosphorylation_type, + phosphorylation_pme)[1]; + conc_test[e, independent_bal_ix] = conc_independent_balanced_experiment; + conc_test[e, dependent_bal_ix] = conc_moiety_pool_test[e]- left_nullspace_independent * conc_independent_balanced_experiment; conc_test[e, unbalanced_mic_ix] = conc_unbalanced_test[e]; vector[N_edge] edge_flux = get_edge_flux(conc_test[e], conc_enzyme_experiment, diff --git a/maud/utility_functions.py b/maud/utility_functions.py index 50d43933d..6f8bb0f82 100644 --- a/maud/utility_functions.py +++ b/maud/utility_functions.py @@ -101,6 +101,54 @@ def get_null_space(a, rtol=1e-5): return v[rank:].T.copy() +def get_left_nullspace(matrix: pd.DataFrame, atol=1e-13, rtol=0.0): + """Compute an approximate basis for the null space (kernel) of a matrix. + + The algorithm used by this function is based on the singular value + decomposition of the given matrix. + + Parameters + ---------- + matrix : ndarray + The matrix should be at most 2-D. A 1-D array with length k + will be treated as a 2-D with shape (1, k) + atol : float + The absolute tolerance for a zero singular value. Singular values + smaller than ``atol`` are considered to be zero. + rtol : float + The relative tolerance for a zero singular value. Singular values less + than the relative tolerance times the largest singular value are + considered to be zero. + + Notes + ----- + If both `atol` and `rtol` are positive, the combined tolerance is the + maximum of the two; that is:: + tol = max(atol, rtol * smax) + Singular values smaller than ``tol`` are considered to be zero. + + Returns + ------- + ndarray + If ``matrix`` is an array with shape (m, k), then the returned + nullspace will be an array with shape ``(k, n)``, where n is the + estimated dimension of the nullspace. + + References + ---------- + Adapted from: + https://scipy.github.io/old-wiki/pages/Cookbook/RankNullspace.html + and then taken from from + https://github.com/opencobra/memote/blob/develop/src/memote/support/consistency_helpers.py#L163 + + """ # noqa: D402 + matrix = np.atleast_2d(matrix) + _, sigma, vh = np.linalg.svd(matrix.T) + tol = max(atol, rtol * sigma[0]) + num_nonzero = (sigma >= tol).sum() + return vh[num_nonzero:].conj().T + + def get_rref(mat): """Return reduced row echelon form of a matrix.""" return sp.Matrix(mat).rref(iszerofunc=lambda x: abs(x) < 1e-10)[0] diff --git a/pyproject.toml b/pyproject.toml index 2496fa957..9429d75cd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ readme = "README.rst" requires-python = ">=3.10" dependencies = [ "pip >= 20", - "arviz >= 0.18.0", + "arviz@git+https://github.com/arviz-devs/arviz#egg=39eddd6", "importlib_resources >= 3.2", "numpy", "scipy", @@ -32,7 +32,7 @@ dependencies = [ "pandas", "matplotlib", "toml", - "cmdstanpy >= 1.2.1", + "cmdstanpy>=1.2.3", "click", "depinfo == 1.7.0", "pydantic >= 2.0", @@ -118,6 +118,8 @@ markers = ["raises"] # Enable the following rules: # pycodestyle (`E`), Pyflakes (`F`), pycodestyle (`W`), flake8-bugbear (`B`) lint.select = ["E", "F", "B", "W", "D"] +# Enable the following rules: +# pycodestyle (`E`), Pyflakes (`F`), pycodestyle (`W`), flake8-bugbear (`B`) lint.ignore = [ "B905", # Use zip() without a `strict` parameter. "D107", # Init methods can be undocumented. @@ -125,7 +127,11 @@ lint.ignore = [ "D213", # https://beta.ruff.rs/docs/rules/multi-line-summary-second-line/ "D104", # __init__.py can be empty. ] +# Enable the following rules: +# pycodestyle (`E`), Pyflakes (`F`), pycodestyle (`W`), flake8-bugbear (`B`) lint.fixable = ["ALL"] +# Enable the following rules: +# pycodestyle (`E`), Pyflakes (`F`), pycodestyle (`W`), flake8-bugbear (`B`) lint.unfixable = [] line-length = 80 target-version = "py310"