diff --git a/maud/getting_stan_inputs.py b/maud/getting_stan_inputs.py index 12da72ea..d1daf2a3 100644 --- a/maud/getting_stan_inputs.py +++ b/maud/getting_stan_inputs.py @@ -135,6 +135,7 @@ def get_network_properties_input( metabolite_codes = codify_maud_object(kinetic_model.metabolites) mic_codes = codify_maud_object(mics) enzyme_codes = codify_maud_object(enzymes) + 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: @@ -187,6 +188,9 @@ def get_network_properties_input( enzyme_codes[er.enzyme_id] if isinstance(er, EnzymeReaction) else 0 for er in edges ] + edge_er_code = [ + er_codes[er.id] if isinstance(er, EnzymeReaction) else 0 for er in edges + ] edge_drain_code = [ drain_codes[d.id] if isinstance(d, Reaction) else 0 for d in edges ] @@ -355,6 +359,7 @@ def get_network_properties_input( "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), "N_pme": len(pme_codes), "N_competitive_inhibition": len(ci_ix_long), @@ -374,6 +379,7 @@ def get_network_properties_input( "ci_mic_ix": ci_mic_codes, "edge_type": edge_mechanism_code, "edge_to_enzyme": edge_enzyme_code, + "edge_to_er": edge_er_code, "edge_to_tc": edge_tc_code, "edge_to_drain": edge_drain_code, "edge_to_reaction": edge_reaction_code, diff --git a/maud/stan/functions.stan b/maud/stan/functions.stan index 830ff928..02540b0f 100644 --- a/maud/stan/functions.stan +++ b/maud/stan/functions.stan @@ -297,12 +297,12 @@ functions { } vector get_vmax_by_edge(vector enzyme, vector kcat, - array[] int edge_to_enzyme, array[] int edge_type) { + array[] int edge_to_enzyme, array[] int edge_to_er, array[] int edge_type) { int N_edge = size(edge_to_enzyme); vector[N_edge] out = rep_vector(1, N_edge); for (f in 1 : N_edge) { if (edge_type[f] != 3) { - out[f] = enzyme[edge_to_enzyme[f]] * kcat[edge_to_enzyme[f]]; + out[f] = enzyme[edge_to_enzyme[f]] * kcat[edge_to_er[f]]; } } return out; @@ -313,7 +313,7 @@ functions { 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 edge_to_enzyme, array[] int edge_to_er, 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, @@ -332,7 +332,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, + vector[N_edge] vmax = get_vmax_by_edge(enzyme, kcat, edge_to_enzyme, edge_to_er, edge_type); vector[N_edge] reversibility = get_reversibility(dgr, temperature, S, conc, edge_type); @@ -371,6 +371,18 @@ functions { sub_by_edge_bounds, edge_type, drain_small_conc_corrector); + // Paste these lines into the function get_edge_flux to debug. + print("free_enzyme_ratio: ", free_enzyme_ratio); + print("vmax: ", vmax); + print("kcat: ", kcat); + print("reversibility: ", reversibility); + print("saturation: ", saturation); + print("allostery: ", allostery); + print("phosphorylation: ", phosphorylation); + print("drain_by_edge: ", drain_by_edge); + print("drain_by_edge: ", drain_by_edge); + print("conc: ", conc); + print("flux: ", vmax .* saturation .* reversibility .* allostery.* phosphorylation .* drain_by_edge); return vmax .* saturation .* reversibility .* allostery .* phosphorylation .* drain_by_edge; } @@ -383,7 +395,7 @@ functions { vector conc_pme, vector drain, real temperature, 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_drain, + array[] int edge_to_enzyme, array[] int edge_to_er, 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, @@ -412,7 +424,7 @@ functions { temperature, drain_small_conc_corrector, S, subunits, edge_type, - edge_to_enzyme, edge_to_drain, + edge_to_enzyme, edge_to_er, edge_to_drain, ci_mic_ix, sub_km_ix_by_edge_long, sub_km_ix_by_edge_bounds, @@ -445,6 +457,7 @@ functions { matrix S, 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, array[] int sub_km_ix_by_edge_long, @@ -469,6 +482,7 @@ functions { int N_edge = cols(S); complex_vector[N_edge] vmax = get_complex_vmax_by_edge(enzyme, kcat, edge_to_enzyme, + edge_to_er, edge_type); vector[N_edge] reversibility = get_reversibility(dgr, temperature, S, conc, edge_type); @@ -512,13 +526,13 @@ functions { } complex_vector get_complex_vmax_by_edge(complex_vector enzyme, vector kcat, - array[] int edge_to_enzyme, + array[] int edge_to_enzyme, array[] int edge_to_er, array[] int edge_type) { int N_edge = size(edge_to_enzyme); complex_vector[N_edge] out = rep_vector(1, N_edge); for (f in 1 : N_edge) { if (edge_type[f] != 3) { - out[f] = enzyme[edge_to_enzyme[f]] * kcat[edge_to_enzyme[f]]; + out[f] = enzyme[edge_to_enzyme[f]] * kcat[edge_to_er[f]]; } } return out; @@ -554,6 +568,7 @@ functions { matrix S, 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, array[] int sub_km_ix_by_edge_long, @@ -576,7 +591,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, + 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, @@ -814,7 +829,7 @@ vector maud_ae_system(vector conc_ind, 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 edge_to_enzyme, array[] into edge_to_er, 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, @@ -858,7 +873,7 @@ vector maud_ae_system(vector conc_ind, temperature, drain_small_conc_corrector, S, subunits, edge_type, - edge_to_enzyme, edge_to_drain, + 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, diff --git a/maud/stan/model.stan b/maud/stan/model.stan index 9454a3a1..84f73fe0 100644 --- a/maud/stan/model.stan +++ b/maud/stan/model.stan @@ -14,6 +14,7 @@ data { int N_prod_km; int N_reaction; int N_enzyme; + int N_er; int N_drain; int N_edge; int N_allostery; @@ -30,6 +31,7 @@ data { array[N_competitive_inhibition] int ci_mic_ix; array[N_edge] int edge_type; // 1 = reversible modular rate law, 2 = drain array[N_edge] int edge_to_enzyme; // 0 if drain + array[N_edge] int edge_to_er; // 0 if drain array[N_edge] int edge_to_tc; // 0 if non-allosteric array[N_edge] int edge_to_drain; // 0 if enzyme array[N_edge] int edge_to_reaction; @@ -220,7 +222,7 @@ transformed parameters { drain_small_conc_corrector, S, left_nullspace_independent, subunits, edge_type, - edge_to_enzyme, edge_to_drain, + 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, @@ -251,7 +253,7 @@ transformed parameters { temperature_train[e], drain_small_conc_corrector, S, subunits, edge_type, - edge_to_enzyme, edge_to_drain, + edge_to_enzyme, edge_to_er, edge_to_drain, ci_mic_ix, sub_km_ix_by_edge_long, sub_km_ix_by_edge_bounds, @@ -407,6 +409,7 @@ generated quantities { subunits, edge_type, edge_to_enzyme, + edge_to_er, edge_to_drain, ci_mic_ix, sub_km_ix_by_edge_long, @@ -450,6 +453,7 @@ generated quantities { subunits, edge_type, edge_to_enzyme, + edge_to_er, edge_to_drain, ci_mic_ix, sub_km_ix_by_edge_long, diff --git a/maud/stan/out_of_sample_model.stan b/maud/stan/out_of_sample_model.stan index 8da0cb74..94890a68 100644 --- a/maud/stan/out_of_sample_model.stan +++ b/maud/stan/out_of_sample_model.stan @@ -14,6 +14,7 @@ data { int N_prod_km; int N_reaction; int N_enzyme; + int N_er; int N_drain; int N_edge; int N_allostery; @@ -29,6 +30,7 @@ data { array[N_competitive_inhibition] int ci_mic_ix; array[N_edge] int edge_type; // 1 = reversible modular rate law, 2 = drain array[N_edge] int edge_to_enzyme; // 0 if drain + array[N_edge] int edge_to_er; // 0 if drain array[N_edge] int edge_to_tc; // 0 if non-allosteric array[N_edge] int edge_to_drain; // 0 if enzyme array[N_edge] int edge_to_reaction; @@ -158,7 +160,7 @@ generated quantities { drain_small_conc_corrector, S, left_nullspace_independent, subunits, edge_type, - edge_to_enzyme, edge_to_drain, + 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, @@ -176,7 +178,7 @@ generated quantities { 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, 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, @@ -188,7 +190,7 @@ generated quantities { temperature_test[e], drain_small_conc_corrector, S, subunits, edge_type, - edge_to_enzyme, edge_to_drain, + edge_to_enzyme, edge_to_er, edge_to_drain, ci_mic_ix, sub_km_ix_by_edge_long, sub_km_ix_by_edge_bounds, diff --git a/pyproject.toml b/pyproject.toml index 7d42aeac..9429d75c 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",