Testing version of bayesianSensitivityTuning

[edkerk]

Description of the new feature:

While rudimentary code was already present, the Bayesian kcat sensitivity tuning was not yet functional.

First introduced in the DLKcat paper (doi:10.1038/s41929-022-00798-z), this is an SMC-ABC (sequential Monte Carlo approximate Bayesian computation) approach to adjust kcat values towards obtaining a model that can more accurately replicate provided experimental data.

The code is now available in fix/bayesianTuning. It can be tested by going through full_ecModel/protocol.m, where the new function is explained from line 283 onwards.

Some elements remain uncertain for now:
1.- What should be used as first standard deviations of kcat values?
2.- If kcatStd = 0 for a particular kcat, will this kcat indeed be kept untouched?
3.- What are good values for the hyperparameters? What are good number of models per (first) generation, number of models to be used to define posterior kcat distributions, etc.
4.- The RMSE does not yet converge as nicely as it did in the DLKcat paper, is that due to the changes made to getrSample (DLKcat version) and updateprior (DLKat version)? I do think that the new code is more robust, as the previous code kept the standard deviation but kcat values in log-scale between generations.

I hereby confirm that:

• The new feature is not already in the main branch of the repository.
• A similar issue does not already exist.

[IVANDOMENZAIN]

kinetic_parameters_summary.txt

Here I provide a txt file providing a summary of my statistical analysis of Kcats reported in BRENDA. The analysis was made on an EC number basis (analysing Kcat distributions per EC number). Only those EC numbers with at least 3 reported Kcats were considered. For each EC number median, mean, maximum, minimum and standard deviation of kcat parameters are calculated.

Additionally, a Shapiro-Wilk test is applied to the log-transformed Kcat distributions per EC numbers, which informs of which Kcat distributions can be fairly considered as log-normal distributed.

Other metrics are also reported, coefficient of variation (CV1) as the ratio between the mean and standard deviation of each distribution, modified coeff. of variation (CV2) as the ratio between the median and standard deviation of each distribution, S = log10(mean)/log10(max/min), which is a measurement of how many orders of magnitude a dsitribution is spanning in relation to its mean value.

Let´s try to use this information to define a proper set of first standard deviations for kcats.

[IVANDOMENZAIN]

One solution could be to assign each kcat with the same standard deviation as that for their corresponding EC number (according to this file), missing EC numbers could be assigned an average std. However this would probably require to upload a new file to the GECKO toolbox (in databases) and retrieve these standard deviation - EC number pairs from it.

[IVANDOMENZAIN]

Another solution, more coarse-grained, would be to assign kcats with standard deviations related to average or global statistic metrics across all reported parameters.

I suggest to asign Kcat with standard deviations as: std = Kcat*/CV (in whixh Kcat* is the prior or initial Kcat and CV is a coefficient of variation representative of the whole kcat dataset.

According to the provided file, the mean CV across all enzyme classes is 1.426, the median is 1.236, and a weighted CV (accounting for the number of entries per enzyme class) is 2.128. Then we could used the highest of these (2.128) as a CV representative of the whole dataset. Therefore, the initial standard deviation for kcats (line 24 in bayesianSensitivityTuning) would be:

kcatStd = ecModel.ec.kcat./2.128;

[IVANDOMENZAIN]

Regarding point number 4, the same for me, RMSE seems to not change across iterations. Let me check the current getrSample more in detail.

[IVANDOMENZAIN]

Related to this, I don´t think RMSE is the most appropiate metric to be used as a distance function in the bayesian approach, as this is not a scaled error metric, which dificults its generalization as a convergence criterion.

RMSE value depends on the magnitude of the variables that are used to compute it. So if very small fluxes are used to compute it (secondary metabolites, for instance), then the resulting RMSE is gonna be a small number, and for high fluxes (carbon source and main byproducts) the RMSE would be high. Even though fluxes are corrected by the carbon content of the exchanged metabolites the resulting RMSE is still not a scaled error metric.

[IVANDOMENZAIN]

As a solution @edkerk and I have discussed that relative changes in RMSE in between iterations can be used as a convergence criterion. For instance, we could add the estimation of relative changes in RMSE as:

RMSE change = (RMSE_t - RMSE_t-1)/RMSE_t-1

Where RMSE_t is the computed RMSE at iteration t and RMSE_t-1 is RMSE at iteration t-1.

Then, this relative RMSE change can be used as a general convergence criterion, as it would not depend on the magnitude of the fluxes used for RMSE estimation.

[IVANDOMENZAIN]

Then a robust convergence criterion could be to stop iterations once the RMSE change is below a given threshold (5% or 1%) for the last n iterations (n around 5, to avoid excessive computation time). This would circumvent the problem of potential noisy behaviour in the asymptotic region of RMSE vs. iteration number.