Model in Wu et al unit test is inconsistent with model in Wu et al notebook

[djinnome]

Hi @djinnome ,

@mcnaughtonadm and @ShantMahserejian met today to work on the BMCA code. Here’s an update of the current status/issue, possible ways forward, and what we will try as an immediate next step:

• The model in the Wu notebook updated to be operation in the test code:

  • Required adding compartments, exchange/sink reactions, and an objective to make the model operational with model.optimize(), which is needed to calculate fluxes
  • The updated model is inconsistent from the original model, especially when considering the Giersch data
    • Newly added reactions/metabolites are missing from the data
    • Prior Trace (and therefore Prior Predictive and FCC) calculations now different
      * Different $E_x$, $E_y$, $N$ matrices
      * Different $v^\ast$
      * Different LinLog model/method needed (vs. original Wu notebook)
  • Error running notebook with missing data for metabs/reactions not included in Giersch data
    * Specifically, stuck on calculating steady state terms: $\chi_{ss}$ and $\hat{v}_ss$.

• Possible alternative approaches:

  1. Continue trying to remedy Wu notebook to work with missing data for metabs/reactions not included in Giersch data
     • Unsure if it will work…
     • Results may be different than what is published, which defeats the purpose of using Wu as a benchmark
  2. Make up numbers for new metabs/reactions by duplicating columns of Giersch data
     • May NOT be ok if FCC-values heavily dependent on other reactions in the model
     • May be ok if FCC-values are isolated to their corresponding reactions $\implies$ just ignore FCCs corresponding to new reactions
       * $\implies$ Likely the case because newly added reaction only handle one metabolite at a time
  3. Use E. coli model in lieu of simple model from Wu paper for testing purposes
     • We are using a method designed around a full genome-scale model, but trying to make it work for a 3-(or 5-)metabolite model
     • Data for E. coli likely to exist, but would be an extra effort at this stage

• Suggest next step – Option 2: Make up numbers
  * We’ve already got a working model for this, so easier step to check if final results agree with original Wu notebook
  * Would be more appropriate rather than changing lines of code impacting the approach
  * Challenge: Unsure of best values to use to augment Gierch data (what values do we use for exchange/sink reactions?)

We’d be happy to chat about this before the weekly PPI meeting if available. We’ll try to keep making progress with the suggest next step to get functional code, and we can revisit how we change the Giersch data after the fact.

[djinnome]

Thanks for the detailed breakdown and analysis.

The updated model is inconsistent from the original model, especially when considering the Giersch data

This is where the work needs to happen. Bayesian MCA uses different conventions than regular FBA, and this pain point has affected the ability of @janisshin to operationalize Bayesian MCA, too.

Recall that there are two different models: A cobra model, and a pymc model. The goal should be to take a cobra model as input and generate a pymc model that exactly reproduces the notebook results. If we can't do this for Wu, then we can't guarantee that it will work for more complex models.

I am happy to help you revise the transformation of the cobra model into a pymc model so that it does reproduce the Wu et al results.

[djinnome]

If this is the input cobra model
https://github.com/pnnl-predictive-phenomics/emll/blob/hackett_methods/tests/test_models/wu2004_model.json

Then this should be the output pymc model:

with pm.Model() as pymc_model:
    
    # Initialize elasticities
    Ex_t = pm.Deterministic('Ex', initialize_elasticity(N, 'ex', b=0.05, sigma=1, alpha=5))
    Ey_t = pm.Deterministic('Ey', initialize_elasticity(-Ey.T, 'ey', b=0.05, sigma=1, alpha=5))
    
    # sample prior
    trace_prior = pm.sample_prior_predictive()
        
    # Error priors. 
    v_err = pm.HalfNormal('v_error', sigma=0.05, initval=.1)
    x_err = pm.HalfNormal('x_error', sigma=0.05, initval=.1)

    # Calculate steady-state concentrations and fluxes from elasticities
    chi_ss, v_hat_ss = ll.steady_state_pytensor(Ex_t, Ey_t, en, yn)

    # Error distributions for observed steady-state concentrations and fluxes
    chi_obs = pm.Normal('chi_obs', mu=chi_ss, sigma=x_err, observed=xn)
    v_hat_obs = pm.Normal('v_hat_obs', mu=v_hat_ss[:, 0].squeeze(),
                          sigma=v_err, observed=vn.squeeze())

    # perform HMC to compute posterior
    trace = pm.sample(500, random_seed=123, tune=1000,step=pm.NUTS(), n_init=1000, nuts_sampler_kwargs=dict(njobs=4))

    # sample posterior
    ppc = pm.sample_posterior_predictive(trace)

Where

• the internal elasticity matrix $E_x$
• the external elasticity matrix $E_y$
• the computed steady-state internal metabolites $\chi_{ss}$
• the computed steady-state fluxes, $\hat{v}_{ss}$
• the observed internal metabolites $\chi_{obs}$
• and the observed fluxes $\hat{v}_{obs}$

are identical to the Wu et al notebook.

[mcnaughtonadm]

So the major problem we need to address here is that our design goal is to try and generalize the workflow of building the pymc model to work for all possible input cobra models. The Wu example is not a standard implementation of this BMCA workflow. It does not utilize a LinLogLeastNorm to obtain the LinLogModel but a LinLogSymbolic2x2 .

import emll
Ex = emll.create_elasticity_matrix(model)

Ey = np.zeros((3, 2))
Ey[0, 0] = 1  # bpg (+) PGM
Ey[2, 1] = 1  # adp (+) PK

ll = emll.LinLogSymbolic2x2(N, Ex, Ey, v_star)

This is a specific behavior only for the Wu "model" (It's not a true cobra model and we have needed to make pretty major changes to it in order to make it into one). This makes it difficult to create a test for BMCA as we would like to use LinLogLeastNorm for all cases as this will be the expected behavior when using a full Genome Scale Model. This is where the hurdle with Wu is. We can create the pymc model no problem, but getting to that point is not the same. If you think we need to dynamically decide which LinLogModel to use, that is an extra layer of complexity that we will have to decide how to handle.

Also this implementation uses a manually defined external elasticity matrix Ey. Which is also something we don't want to be doing as it is specific to Wu and we would not be able to automatically assign it like we do with genome scale models by using create_Ey_matrix . It does seem like Peter St. John's documentation for this method recommends we handle Ey creation manually but I've found it to be pretty accurate for well-defined cobra models, which Wu is not, hence the manually definition.

So we diverge in a few areas:

1. Our workflow is really set up to take in genome scale cobra models, and not 3 reaction incomplete models
2. Wu uses some hand wavy assignments due to it not being a true cobra model which is incompatible with how we restructured the code (defining Ey as the main example).
3. To make a true test case for Wu, we would be using methods outside of our typical use case so we aren't truly testing our new generalized code, but old methods we probably won't use in the end product.

Thoughts after writing all this:

I think we'd just need to step back and look at it again. We could pre-define Ex, Ey for the ll (linlog model) but I feel like that doesn't really help us. Our problem is just that Wu requires a lot of very specific assignments that don't mesh well with the approach we want to take with the workflow, I know it's a "simple" model, but I think we are having the opposite problem that we do with hackett where it's TOO simple and doesn't help actually test our code.

[djinnome]

Hey @mcnaughtonadm I agree that the Wu model contains lots of corner cases, so perhaps from a strategic perspective, it would make sense to punt on Wu for the time being, and focus on more typical workflows such as Contador or Hackett or A. niger so that at least you have one or two examples working end-to-end. We can always return to Wu later once the pipeline has been hardened so that there are fewer moving parts to deal with and we can just focus on what is different between Wu and the typical use case.

[djinnome]

In fact, maybe A. niger would be the best one to work on, since it works in our hands, it has been experimentally validated, and is an example of a "typical" workflow that we would expect to get.

[ShantMahserejian]

@djinnome and @mcnaughtonadm: please see the latest notebook file that I just pushed to the hacket_methods branch since your last comments. I updated the notebook to include a comparison of FCCs sampled from the prior predictive trace between those generated from the original Wu notebook and those generated using the updated cobra model (see violin plots output from cells 18 and 19).

I left in the original code snippets as commented lines where they were relevant. Some of the notebook cell ordering has change slightly as we were exploring things. (Cells 3 and 4 for the cobra model; cell 8 for the elasticity matrices; cell 10 for the Giersch data loading, though the new code here is commented)

---

Short version of findings:
A slight difference in the FCCs sampled from the prior predictive distributions are observed in the violin plots comparing results from the different models and corresponding differences in elasticity matrices. This indicates a departure from the published work even at this early stage of the analysis. Since the changes to the model can be deemed necessary for making progress with the code being developed/tested, this provides further support selecting a different model for establishing a baseline for testing and comparison.

---

Longer version, including recapitulation of the discussion so far:
@djinnome : yes, your code for how the pymc model should be produced is valid, but as the notebook breaks it down into smaller steps, we can see a departure between this generalized approach to the highly specified approach used in the original Wu notebook. The list of issues at the top of this thread were discovered during our tests with generating prior predictive distributions, but we anticipate additional issues with posterior predictives as well. Please note the following breakdown of these ideas:

• The cobra model is updated from the original notebook version to now include compartments, exchange/sink reactions, and an objective which are all needed to run model.optimize(). This is needed to calculate fluxes and satisfy SBML requirements for saving the file, without which we cannot test FCC calculations using the "Wu" model.

• Using only the cobra model (and without having to run model.optimize()), we can calculate the internal and external elasticity matrices ($E_x$ and $E_y$) and the stoichiometric matrix, $N$. This is enough generate an unconditional pymc model, which can be used to sample a prior predictive trace, trace_prior. (Please see build_unconditional_pymc_model() in the pymc_utils.py in our code, and cell # 8 in the notebook for how we 've had to modify these calculations with the new model, which are effectively the same as the code shared above). ==> Problem # 1: the $E_x$ and $E_y$ matrices already look different from those corresponding to the original cobra model, so we can expect trace_prior to also be impacted. Namely the matrix shapes are different because of the added metabolites and reactions in the cobra model

• Next, we also need an array of fluxes, v_star which we now obtain by running model.optimize(), whereas the original Wu notebook uses the reference strain flux from the Giersch data to calculate their version of v_star ==> Problem # 2: to produce a truly unconditional prior predictive distribution (i.e. a bare minimum approach that is blind to any data), we need to change how the Wu notebook calculates v_star, at least as an initial step to test our more complete pipeline (see cell # 7 in the notebook, which is identical to how util.calculate_v_star() calculates v_star in sample_prior_fccs() in the pymc_utils.py file).

• In order to calculate FCCs from the prior predictive distribution, we need to generate a LinLog model using $E_x$, $E_y$, $N$ and v_star. ==> Problem # 3: due to the differences in the cobra model and the corresponding calculated matrices/arrays, the LinLog model needs to be created using a different method than the original Wu notebook (the matrix shapes are different, so emll.LinLogSymbolic2x2()) (see how ll is defined in cell # 8, which is now similar to how it is defined in the sample_prior_fccs() function in pymc_utils.py).

• Now that the FCCs can be sampled the prior predictives, we can compare the results between the original Wu notebook and those resulting from using a modified cobra model and downstream calculations. ==> Problem # 4: The different versions of FCC corresponding to the 3 reactions of interest (not exchange/sink reactions that we were forced to add) produce plots with slight differences in centers and variance (see plots generated from cells # 18 and # 19). it's difficult to attribute how much of the differences are due to random sampling and how much is due to differences in cobra models, corresponding matrices/arrays, and LinLog models.

• In order to continue with sampling the posterior distribution conditioned on the Giersch data, we need to first construct a probabilistic model, but we encounter an error when trying to calculate the steady state terms chi_ss and v_hat_ss (see cell # 20). ==> Problem # 5: we cannot move forward using the Giersch data as-is, since it is incompatible with the updated version of the cobra model and and the corresponding differences in shapes of the elasticity matrices. Modifying the Giersch data to accommodate the changes due to additional metabolites and reactions in the cobra model is even further departure from the published Wu results, and making it less of a viable option for use as a baseline case for testing and comparison. (see commented lines in cell # 10, which are possible logical changes to the Giersch data that can be used, but they have not been validated to remove the current errors seen with cell # 20).