Skip to content

Commit

Permalink
fix: introduce right indexing for kcats (edge to enz reaction) (#465)
Browse files Browse the repository at this point in the history 
* feat: left_nullspace for conserved moieties

* fix: introduce a edge_to_er mapping for kcats

* fix: remove Nick's lingering left nullspace

* test: update test expected inputs

* refactor: use proper plurals for er indices

* style: ruff format

---------

Co-authored-by: NicholasCowie <[email protected]>
  • Loading branch information
@carrascomj @NicholasCowie
carrascomj and NicholasCowie authored Jun 21, 2024
1 parent 99ecc35 commit 9747fe6
Show file tree
Hide file tree
Showing 6 changed files with 37 additions and 14 deletions.
2 changes: 1 addition & 1 deletion maud/data/example_inputs/example_ode/input_data_train.json
View file
Original file line number Diff line number Diff line change
@@ -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": [], "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_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]]}
6 changes: 6 additions & 0 deletions maud/data/example_inputs/linear/expected_stan_input.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
"N_edge": 3,
"N_unbalanced": 2,
"N_enzyme": 3,
"N_er": 3,
"N_phosphorylation": 0,
"N_pme": 0,
"N_competitive_inhibition": 1,
Expand Down Expand Up @@ -59,6 +60,11 @@
2,
3
],
"edge_to_er": [
1,
2,
3
],
"edge_to_tc": [
1,
2,
Expand Down
6 changes: 6 additions & 0 deletions maud/getting_stan_inputs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,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 allosteries is not None:
Expand Down Expand Up @@ -177,6 +178,9 @@ def get_network_properties_input(
enzyme_codes[er.enzyme_id] if isinstance(er, EnzymeReaction) else 0
for er in edges
]
edge_er_codes = [
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
]
Expand Down Expand Up @@ -338,6 +342,7 @@ def get_network_properties_input(
"N_edge": len(S.columns),
"N_unbalanced": len(unbalanced_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),
Expand All @@ -355,6 +360,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_codes,
"edge_to_tc": edge_tc_code,
"edge_to_drain": edge_drain_code,
"edge_to_reaction": edge_reaction_code,
Expand Down
23 changes: 14 additions & 9 deletions maud/stan/functions.stan
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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,
Expand All @@ -333,6 +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);
vector[N_edge] reversibility = get_reversibility(dgr, temperature, S,
conc, edge_type);
Expand Down Expand Up @@ -382,7 +383,7 @@ functions {
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,
Expand Down Expand Up @@ -410,7 +411,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,
Expand Down Expand Up @@ -443,6 +444,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,
Expand All @@ -467,6 +469,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);
Expand Down Expand Up @@ -511,12 +514,13 @@ functions {

complex_vector get_complex_vmax_by_edge(complex_vector enzyme, vector kcat,
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;
Expand Down Expand Up @@ -552,6 +556,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,
Expand All @@ -574,7 +579,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,
Expand Down Expand Up @@ -810,7 +815,7 @@ vector maud_ae_system(vector conc,
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,
Expand Down Expand Up @@ -843,7 +848,7 @@ vector maud_ae_system(vector conc,
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,
Expand Down
8 changes: 6 additions & 2 deletions maud/stan/model.stan
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ data {
int<lower=1> N_prod_km;
int<lower=1> N_reaction;
int<lower=1> N_enzyme;
int<lower=1> N_er;
int<lower=0> N_drain;
int<lower=1> N_edge;
int<lower=0> N_allostery;
Expand All @@ -25,6 +26,7 @@ data {
array[N_competitive_inhibition] int<lower=1, upper=N_mic> ci_mic_ix;
array[N_edge] int<lower=1, upper=3> edge_type; // 1 = reversible modular rate law, 2 = drain
array[N_edge] int<lower=0, upper=N_enzyme> edge_to_enzyme; // 0 if drain
array[N_edge] int<lower=0, upper=N_er> edge_to_er; // 0 if drain
array[N_edge] int<lower=0, upper=N_allostery> edge_to_tc; // 0 if non-allosteric
array[N_edge] int<lower=0, upper=N_drain> edge_to_drain; // 0 if enzyme
array[N_edge] int<lower=0, upper=N_reaction> edge_to_reaction;
Expand Down Expand Up @@ -201,7 +203,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,
prod_km_ix_by_edge_long,
Expand Down Expand Up @@ -231,7 +233,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,
Expand Down Expand Up @@ -388,6 +390,7 @@ generated quantities {
subunits,
edge_type,
edge_to_enzyme,
edge_to_er,
edge_to_drain,
ci_mic_ix,
sub_km_ix_by_edge_long,
Expand Down Expand Up @@ -431,6 +434,7 @@ generated quantities {
subunits,
edge_type,
edge_to_enzyme,
edge_to_er,
edge_to_drain,
ci_mic_ix,
sub_km_ix_by_edge_long,
Expand Down
Loading

0 comments on commit 9747fe6

Please sign in to comment.