diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 879e275a..b594fd61 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -136,3 +136,14 @@ @article{kacser1973 langid = {english}, keywords = {Arginine,Enzyme Repression,Enzymes,Feedback,Kinetics,Melanins,Metabolism,{Models, Chemical},Neurospora crassa,Tyrosine 3-Monooxygenase} } +@inproceedings{margossianComputingSteadyStates2018, + title = {Computing Steady States with {{Stan}}'s Nonlinear Algebraic Solver}, + booktitle = {Stan {{Conference}} 2018}, + author = {Margossian, Charles}, + year = {2018}, + month = jan, + doi = {10.5281/zenodo.1284375}, + url = {10.5281/zenodo.1284375}, + abstract = {Stan's numerical algebraic solver can be used to solve systems of nonlinear algebraic equations with no closed form solutions. One of its key applications in scientific and engineering fields is the computation of equilibrium states (equivalently steady states). This case study illustrates the use of the algebraic solver by applying it to a problem in pharmacometrics. In particular, I show the algebraic system we solve can be quite complex and embed, for instance, numerical solutions to ordinary differential equations. The code in R and Stan are provided, and a Bayesian model is fitted to simulated data}, + file = {/Users/tedgro/Zotero/storage/5XCZ7KVA/Margossian - 2018 - Computing steady states with Stan's nonlinear alge.pdf} +} diff --git a/docs/inputting.md b/docs/inputting.md index c0b3f820..b15464f9 100644 --- a/docs/inputting.md +++ b/docs/inputting.md @@ -33,7 +33,8 @@ The following optional fields can also be specified: * `reject_non_steady` Boolean saying whether to reject draws that enter non-steady states * `penalize_non_steady` Boolean saying whether to penalize steady state deviations in the likelihood. It cannot be `True` when `reject_non_steady` is `True`. -* `ode_config` Table of configuration options for Stan's ode solver +* `ode_solver_config` Table of configuration options for Stan's ode solver +* `algebra_solver_config` Table of configuration options for Stan's algebra solver * `cmdstanpy_config` Table of keyword arguments to the cmdstanpy method [`CmdStanModel.sample](https://cmdstanpy.readthedocs.io/en/v1.1.0/api.html#cmdstanpy.CmdStanModel.sample) * `cmdstanpy_config_predict` Table of overriding sample keyword argments for predictions * `stanc_options` Table of valid choices for [CmdStanModel](https://cmdstanpy.readthedocs.io/en/v1.1.0/api.html#cmdstanpy.CmdStanModel) argument `stanc_options` @@ -68,11 +69,15 @@ Here is an example configuration file: chains = 4 save_warmup = true - [ode_config] + [ode_solver_config] + abs_tol = 1e-4 + rel_tol = 1e-4 + max_num_steps = 1e6 + + [algebra_solver_config] abs_tol = 1e-4 rel_tol = 1e-4 max_num_steps = 1e6 - timepoint = 1e3 This file tells Maud that a file representing a kinetic model can be found at the relative path `kinetic_model.toml`, and that priors and experimental diff --git a/docs/theory.md b/docs/theory.md index bb2c1568..8d151963 100644 --- a/docs/theory.md +++ b/docs/theory.md @@ -226,9 +226,8 @@ solving [](#eq-steady). :::{note} -Note that in practice we do not solve [](#eq-steady) directly but instead use -ODE simulation - see section {ref}`sec-solving-the-steady-state-problem` for -details +See section {ref}`sec-solving-the-steady-state-problem` for how we solve +[](#eq-steady) in practice. ::: @@ -337,16 +336,19 @@ adaptive Hamiltonian Monte Carlo requires gradients of the posterior distribution, it is also necessary to calculate sensitivities of the steady state solution with respect to all parameters. -Our approach to this problem is to choose a starting concentration vector -$x_{0}$ and a simulation time $t$, then find $x_t$ using numerical ODE -integration. To verify whether $x_t$ is a steady state we then evaluate $S\cdot -f_v(x_t, \theta)$ and check if the result is sufficiently close to zero. - -Stan provides an interface to the Sundials ODE solver CVODES, including -gradient calculations. - -For the systems we have investigated, this method works better than solving the -steady state problem using an algebra solver. +Our approach to this problem is, before MCMC sampling, to choose a +concentration vector $x_{guess}$. Every time Maud's sampler evaluates the +statistical model's log probability and gradients, it numerically solves +[](#eq-steady), and finds its gradients, using this starting guess and Stan's +interface to the Sundials Newton solver IDAS. + +The system function that is passed to IDAS includes a representation of +[](#eq-steady) in Stan code. In order to model the instantaneous change of +metabolite concentrations for a given concentration vector due to +[](#eq-steady), Maud uses the Stan interface to the Sundials ODE solver CVODES +to evolve the system. This hybrid approach, using both numerical ODE solving +and algebra solving, follows the one taken in +{cite}`margossianComputingSteadyStates2018`. It is relevant to note a caveat for the Maximum A Posteriori estimation (optimization) in opposition to sampling the full distribution. In an diff --git a/maud/data/example_inputs/example_ode/config.toml b/maud/data/example_inputs/example_ode/config.toml index 2d828ba2..1178482a 100644 --- a/maud/data/example_inputs/example_ode/config.toml +++ b/maud/data/example_inputs/example_ode/config.toml @@ -9,6 +9,5 @@ iter_warmup = 0 iter_sampling = 1 chains = 1 -[ode_config] +[ode_solver_config] max_num_steps = 1e9 -timepoint = 1e3 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 bbf0ffae..ea6cbd44 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_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_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": [], "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_drain_test": [[[]], [[]]], "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_enzyme_test": [[[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_unbalanced_test": [[[1.6094379124341003, -0.6931471805599453]], [[0.2, 0.2]]], "priors_conc_pme_train": [[[]], [[]]], "priors_conc_pme_test": [[[]], [[]]], "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": 1e-09, "abs_tol": 1e-09, "timepoint": 1000.0, "max_num_steps": 1000000000, "conc_init": [[0.323117, 3.02187]]} +{"N_mic": 4, "N_edge_sub": 5, "N_edge_prod": 5, "N_edge": 5, "N_unbalanced": 2, "N_enzyme": 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_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": [], "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]]} diff --git a/maud/data/example_inputs/linear/config.toml b/maud/data/example_inputs/linear/config.toml index f2f8414d..47dd698c 100644 --- a/maud/data/example_inputs/linear/config.toml +++ b/maud/data/example_inputs/linear/config.toml @@ -13,8 +13,12 @@ chains = 4 save_warmup = true seed = 1234 -[ode_config] +[ode_solver_config] +abs_tol = 1e-7 +rel_tol = 1e-7 +max_num_steps = 1e6 + +[algebra_solver_config] abs_tol = 1e-5 rel_tol = 1e-5 max_num_steps = 1e6 -timepoint = 1e3 diff --git a/maud/data/example_inputs/linear/expected_stan_input.json b/maud/data/example_inputs/linear/expected_stan_input.json index 889d005b..f08ca74e 100644 --- a/maud/data/example_inputs/linear/expected_stan_input.json +++ b/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": 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, "penalize_non_steady": 0, "steady_state_threshold_abs": 1e-06, "steady_state_penalty_rel": 1e-08, "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, "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-07, "abs_tol_ode": 1e-07, "max_num_steps_ode": 1000000, "rel_tol_alg": 1e-05, "abs_tol_alg": 1e-05, "max_num_steps_alg": 1000000, "conc_init": [[0.59, 0.38], [0.54, 0.38]]} diff --git a/maud/data/example_inputs/methionine/config.toml b/maud/data/example_inputs/methionine/config.toml index a745cb6c..7ec0eda4 100644 --- a/maud/data/example_inputs/methionine/config.toml +++ b/maud/data/example_inputs/methionine/config.toml @@ -9,15 +9,19 @@ reject_non_steady = true [cmdstanpy_config] iter_warmup = 1000 iter_sampling = 1000 -max_treedepth = 11 +max_treedepth = 10 chains = 4 save_warmup = true refresh = 1 metric = "dense_e" adapt_delta = 0.99 -[ode_config] -rel_tol = 1e-9 -abs_tol = 1e-9 +[ode_solver_config] +rel_tol = 1e-8 +abs_tol = 1e-8 max_num_steps = 1000000.0 -timepoint = 1e+20 + +[algebra_solver_config] +rel_tol = 1e-6 +abs_tol = 1e-6 +max_num_steps = 1000000 diff --git a/maud/data_model/maud_config.py b/maud/data_model/maud_config.py index 22707a02..dfe300de 100644 --- a/maud/data_model/maud_config.py +++ b/maud/data_model/maud_config.py @@ -1,18 +1,27 @@ """Provides dataclass MaudConfig.""" from typing import Optional +from pydantic import Field from pydantic.class_validators import root_validator -from pydantic.dataclasses import Field, dataclass +from pydantic.dataclasses import dataclass @dataclass(frozen=True) -class ODEConfig: - """Config that is specific to the ODE solver.""" +class ODESolverConfig: + """Config that is specific to an ODE solver.""" rel_tol: float = 1e-9 abs_tol: float = 1e-9 - max_num_steps: int = int(1e9) - timepoint: float = 500 + max_num_steps: int = int(1e7) + + +@dataclass(frozen=True) +class AlgebraSolverConfig: + """Config that is specific to an ODE solver.""" + + rel_tol: float = 1e-7 + abs_tol: float = 1e-7 + max_num_steps: int = int(1e6) @dataclass @@ -27,7 +36,8 @@ class MaudConfig: :param cmdstanpy_config: Arguments to cmdstanpy.CmdStanModel.sample. :param reject_non_steady: Reject draws if a non-steady state is encountered. :param penalize_non_steady: Penalize the deviation from steady state in the log likelihood. - :param ode_config: Configuration for Stan's ode solver. + :param ode_solver_config: Configuration for Stan's ode solver. + :param algebra_solver_config: Configuration for Stan's algebra solver. :param stanc_options: Options for CmdStanModel argument `stanc_options`. :param cpp_options: Options for CmdStanModel `cpp_options`. :param variational_options: Arguments for CmdStanModel.variational. @@ -55,7 +65,10 @@ class MaudConfig: variational_options: Optional[dict] = None optimize_options: Optional[dict] = None user_inits_file: Optional[str] = None - ode_config: ODEConfig = Field(default_factory=ODEConfig) + ode_solver_config: ODESolverConfig = Field(default_factory=ODESolverConfig) + algebra_solver_config: AlgebraSolverConfig = Field( + default_factory=AlgebraSolverConfig + ) reject_non_steady: bool = True penalize_non_steady: bool = False steady_state_threshold_abs: float = 1e-8 diff --git a/maud/getting_stan_inputs.py b/maud/getting_stan_inputs.py index ce450d23..52eddef6 100644 --- a/maud/getting_stan_inputs.py +++ b/maud/getting_stan_inputs.py @@ -467,10 +467,12 @@ def get_config_input(config: MaudConfig): "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, - "rel_tol": config.ode_config.rel_tol, - "abs_tol": config.ode_config.abs_tol, - "timepoint": config.ode_config.timepoint, - "max_num_steps": config.ode_config.max_num_steps, + "rel_tol_ode": config.ode_solver_config.rel_tol, + "abs_tol_ode": config.ode_solver_config.abs_tol, + "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, + "max_num_steps_alg": config.algebra_solver_config.max_num_steps, } diff --git a/maud/stan/functions.stan b/maud/stan/functions.stan index 90d2320c..016ba938 100644 --- a/maud/stan/functions.stan +++ b/maud/stan/functions.stan @@ -793,4 +793,66 @@ functions { } return out; } + +vector maud_ae_system(vector conc, + data real rel_tol_ode, data real abs_tol_ode, data int max_num_steps_ode, + vector conc_unbalanced, + 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, + vector conc_pme, vector drain, real temperature, + real drain_small_conc_corrector, matrix S, + vector subunits, array[] int edge_type, + array[] int edge_to_enzyme, array[] int edge_to_drain, + array[] int ci_mic_ix, + array[] int sub_km_ix_by_edge_long, + array[,] int sub_km_ix_by_edge_bounds, + array[] int prod_km_ix_by_edge_long, + array[,] int prod_km_ix_by_edge_bounds, + array[] int sub_by_edge_long, + array[,] int sub_by_edge_bounds, + array[] int prod_by_edge_long, + array[,] int prod_by_edge_bounds, + array[] int ci_ix_long, array[,] int ci_ix_bounds, + array[] int allostery_ix_long, + array[,] int allostery_ix_bounds, + array[] int allostery_type, array[] int allostery_mic, + array[] int edge_to_tc, + array[] int phosphorylation_ix_long, + 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}, + rel_tol_ode, abs_tol_ode, max_num_steps_ode, + conc_unbalanced, + balanced_mic_ix, + unbalanced_mic_ix, + conc_enzyme, + dgr, kcat, km, ki, + tc, + dc, kcat_pme, + conc_pme, + drain, + temperature, + drain_small_conc_corrector, S, + subunits, edge_type, + edge_to_enzyme, 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]); +} + } diff --git a/maud/stan/model.stan b/maud/stan/model.stan index 02ead00b..9ca60e4e 100644 --- a/maud/stan/model.stan +++ b/maud/stan/model.stan @@ -90,15 +90,17 @@ data { array[2, N_experiment_train] vector[N_drain] priors_drain_train; // configuration array[N_experiment_train] vector[N_mic - N_unbalanced] conc_init; - real rel_tol; - real abs_tol; + 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 max_num_steps; int likelihood; // set to 0 for priors-only mode real drain_small_conc_corrector; - real timepoint; int reject_non_steady; int penalize_non_steady; } @@ -154,7 +156,7 @@ transformed parameters { flux_train[e] = rep_vector(0, N_reaction); vector[N_enzyme] conc_enzyme_experiment = conc_enzyme_train[e]; vector[N_pme] conc_pme_experiment = conc_pme_train[e]; - array[1] vector[N_mic - N_unbalanced] conc_balanced_experiment; + vector[N_mic - N_unbalanced] conc_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) { @@ -170,39 +172,42 @@ transformed parameters { e); conc_pme_experiment[pko_experiment] = rep_vector(0, N_pko_experiment); } - conc_balanced_experiment = ode_bdf_tol(dbalanced_dt, conc_init[e], - initial_time, {timepoint}, - rel_tol, abs_tol, max_num_steps, - 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_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[1]; + 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_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]; { vector[N_edge] edge_flux = get_edge_flux(conc_train[e], @@ -238,42 +243,6 @@ transformed parameters { for (j in 1 : N_edge) { flux_train[e, edge_to_reaction[j]] += edge_flux[j]; } - if (reject_non_steady == 1 - && check_steady_state((S * edge_flux)[balanced_mic_ix], - conc_balanced_experiment[1], - steady_state_threshold_abs, - steady_state_threshold_rel) - == 0) { - print("Non-steady state in experiment ", e); - print("Balanced metabolite concentration", - conc_balanced_experiment[1]); - print("flux_train: ", flux_train[e]); - print("conc_init: ", conc_init); - print("conc_unbalanced_train: ", conc_unbalanced_train[e]); - print("log_conc_unbalanced_train_z: ", log_conc_unbalanced_train_z[e]); - print("conc_enzyme_experiment: ", conc_enzyme_experiment); - print("log_conc_enzyme_train_z: ", log_conc_enzyme_train_z[e]); - print("km: ", km); - print("log_km_z: ", log_km_z); - print("drain_train: ", drain_train[e]); - print("drain_train_z: ", drain_train_z[e]); - print("kcat: ", kcat); - print("log_kcat_z: ", log_kcat_z); - print("dgr_train: ", dgr_train[e]); - print("ki: ", ki); - print("log_ki_z: ", log_ki_z); - print("dissociation_constant: ", dissociation_constant); - print("log_dissociation_constant_z: ", log_dissociation_constant_z); - print("transfer_constant: ", transfer_constant); - print("log_transfer_constant_z: ", log_transfer_constant_z); - print("kcat_pme: ", kcat_pme); - print("log_kcat_pme_z: ", log_kcat_pme_z); - print("conc_pme_experiment: ", conc_pme_experiment); - print("log_conc_pme_train_z: ", log_conc_pme_train_z[e]); - print("psi_train: ", psi_train); - print("psi_train_z: ", psi_train_z); - reject("Rejecting"); - } } } }