kcat sensitivity tuning#

MATLAB functions in GECKO/src/kcat_tuning of the GECKO toolbox.

Functions#

Function	Summary
`bestParsimonious`	The fewest-changes point whose RMSE is within
`cmaesKcatTuning`	Tune kcats against experimental data with CMA-ES.
`evotuneScore`	Score one ecModel's current kcats against evotune data.
`fillCarbonNum`	Add a carbon number field to the model.
`isozymeTieMap`	Index of the parameter each isozyme copy follows.
`loadEvotuneData`	Load experimental data for evotune (CMA-ES) kcat tuning.
`nChanged`	How many kcats differ from their prior by more than relTol.
`parsimonyFrontier`	Score kcat with progressively more of it reverted to
`resolveTunableContext`	The tunable kcat subset plus everything derived
`screenKcatLeverage`	Report which kcats the data can actually speak to.
`selectTunableMask`	Build a tunableMask from a screenKcatLeverage report.
`sensitivityTuning`	Relax the most limiting kcats to reach a growth rate.
`sigmaFitter`	Fit the average enzyme saturation factor of an ecModel.
`truncateValues`	Round values in selected columns to six significant digits.
`tunePriorPenaltyWeight`	Sweep priorPenaltyWeight and report the

Reference#

bestParsimonious#

function point = bestParsimonious(points,varargin)

The fewest-changes point whose RMSE is within

tolerance of the best on a parsimonyFrontier.

Defaults to 2%, which is the scale of alternate-optimum noise typically seen from re-solving an LP with a slightly different kcat vector -- below that, two vectors are not meaningfully different in fit, so the sparser one wins.

Parameters

points : struct array: as returned by parsimonyFrontier.

Name-Value Arguments

tolerance : double: default 0.02.

Returns

point : struct: the chosen entry of points.

See also

parsimonyFrontier

cmaesKcatTuning#

function [ecModel,result] = cmaesKcatTuning(ecModel,varargin)

Tune kcats against experimental data with CMA-ES.

Configuration comes entirely from modelAdapter.params.evotune: trust tiers (kcatSources/sigma0logSource/sigma0logDefault), maxGrowthWeight, priorPenaltyWeight and tieIsozymes control the search, and maxGenerations/rmseThreshold are its stopping conditions.

Vendors a separable (diagonal-covariance) CMA-ES -- MATLAB has no built-in implementation and Global Optimization Toolbox doesn't ship one either. Separable rather than full-covariance because that is what geckopy's own port actually runs (cma's CMA_diagonal=True) and what its defaults were validated against; it is also lower risk to implement correctly. See the "CMA-ES core" local functions below.

Parameters

ecModel : struct: ecModel with a populated ec.kcat.

Name-Value Arguments

modelAdapter : ModelAdapter: a loaded model adapter (default: the current default model adapter).
evoData : struct: as returned by loadEvotuneData (default: loaded from modelAdapter).
tunableMask, screen, targetImpactShare: control which kcats are searched over, in order of precedence. Pass tunableMask to fix the set yourself. Otherwise, a mask is built by selectTunableMask at targetImpactShare (default 0.9), from screen if given, or from a fresh screenKcatLeverage call otherwise.
popsize : double: CMA-ES population size. Defaults to 4 + floor(3*log(n)), the standard dimension-scaled default, so it adapts to however many free parameters the mask/tying produced.
nProc : double: number of parallel workers used to score each generation's candidates (default 1, i.e. serial). Requires Parallel Computing Toolbox for > 1.
seed : double: if given, seeds the random draws for a reproducible run.
verbose : logical: whether per-generation progress is printed (default true).

Returns

ecModel : struct: mutated: tunable rows carry the best kcat vector CMA-ES found, with kcat constraints already applied. MATLAB structs are copied on assignment, so -- unlike geckopy's in-place mutation -- this output must be captured by the caller.
result : struct: rxns/oldKcat/newKcat/groups cover the selected tunable set, tied groups included (their members share one value in newKcat). rmseTrace/objectiveTrace record the best-so-far plain RMSE and optimised objective per generation. nGenerations is how many generations actually ran. converged is true if rmseTrace(end) <= params.rmseThreshold.

Raises

error: if there are no tunable kcats, or fewer than two free parameters remain after masking and tying -- too little for CMA-ES to search over.

See also

screenKcatLeverage, selectTunableMask, tunePriorPenaltyWeight

evotuneScore#

function [objective,rmse,rmseDetail,ecModel] = evotuneScore(ecModel,evoData,varargin)

Score one ecModel's current kcats against evotune data.

Internal helper used by cmaesKcatTuning, screenKcatLeverage and tunePriorPenaltyWeight to score one ecModel (with a given set of kcats already applied via applyKcatConstraints) against the experimental growth and flux data. Fuses what geckopy splits into simulate.py (the FBA half) and distance.py (the RMSE half) into one function, matching this module's own predecessor abc_max.m -- there is no MATLAB test suite here to benefit from Python's testability-driven split.

Parameters

ecModel : struct: ecModel with the current candidate kcats already applied (via applyKcatConstraints). Not mutated.
evoData : struct: structure with experimental data to be used, as returned by loadEvotuneData.

Name-Value Arguments

modelAdapter : ModelAdapter: a loaded model adapter (default: the current default model adapter).
maxGrowthWeight : double: relative weight of the max-growth dataset against the flux dataset: (rmseFlux + w*rmseMaxGrate) / (w+1) (default 1.0).
priorPenaltyWeight : double: weight on a prior term added to the objective (default 0, i.e. the objective equals the plain RMSE). See kcatVec/kcat0/sigma0log below.
kcatVec, kcat0, sigma0log : double: only used when priorPenaltyWeight > 0: the tunable kcat subset's current values, prior values, and per-kcat prior log-space std dev, all the same length and in the same order. The penalty added to the objective is priorPenaltyWeight * mean((log(kcatVec./kcat0)./sigma0log).^2).

Returns

objective : double: rmse, plus the prior penalty term when priorPenaltyWeight > 0. This is what a search should minimise.
rmse : double: the plain (unpenalised) average RMSE, always comparable across runs regardless of priorPenaltyWeight.
rmseDetail : cell: per-condition RMSE values, labelled by data type and condition index.
ecModel : struct: the input ecModel, with ec.excarbon cached on it if it wasn't already present. MATLAB structs are copied on assignment, not shared by reference like a Python object, so callers scoring many candidates in a loop should carry this output back in as the base copy for the next candidate, rather than recomputing excarbon (an expensive metabolite-formula scan) on every single call.

See also

cmaesKcatTuning, screenKcatLeverage, loadEvotuneData

fillCarbonNum#

function model = fillCarbonNum(model)

Add a carbon number field to the model.

Adds a model.excarbon field with the carbon number of each exchange reaction, used to weight the RMSE during evotune kcat tuning.

Parameters

model : struct: ecModel that was generated by makeEcModel, or loaded from an earlier run. Not compatible with ecModels generated by earlier GECKO versions (pre 3.0).

Returns

model : struct: ecModel with an added model.excarbon field.

isozymeTieMap#

function tieMap = isozymeTieMap(rxnIds,kcat0,varargin)

Index of the parameter each isozyme copy follows.

A reaction catalysed by several isozymes gets one ec.rxns row per isozyme, each with its own kcat. Where the kcat assignment could tell them apart, those are genuinely different parameters; where it could not, every copy shares the same value from the same source, and an untied search is then free to invent a difference between them that means nothing biologically. Copies group together (and share one parameter) when they share a base reaction (the part of the rxn id before its "EXP" isozyme suffix), a prior value, and -- when sources is given -- a source.

Parameters

rxnIds : cell: ec.rxns entries.
kcat0 : double: prior kcat per row, same length as rxnIds.

Name-Value Arguments

sources : cell: ec.source entries, same length as rxnIds (default: not used, so grouping ignores source).
relTol : double: relative tolerance for treating two prior values as equal (default 1e-9).

Returns

tieMap : double: one entry per position. tieMap(i) == i for a free parameter, otherwise the (lower) index of the representative it follows.

See also

cmaesKcatTuning, screenKcatLeverage

loadEvotuneData#

function evoData = loadEvotuneData(varargin)

Load experimental data for evotune (CMA-ES) kcat tuning.

Reads the three experimental data files used by cmaesKcatTuning from the data subfolder of the model directory defined by the model adapter. Each file is optional: a missing file yields [] (fluxData/maxGrate) or an empty cell array (zeroFlux) rather than an error.

Name-Value Arguments

modelAdapter : ModelAdapter: a loaded model adapter (default: the current default model adapter).

Returns

evoData : struct: structure with the loaded experimental data.

Notes

The evoData structure has the following fields:

fluxData : flux data loaded from evotuneFluxData.tsv, with the biomass reaction id stored in fluxData.biomass. [] if that file isn't present.
maxGrate : maximum growth-rate data loaded from evotuneMaxGrowth.tsv, with the biomass reaction id stored in maxGrate.biomass. [] if that file isn't present.
zeroFlux : reaction IDs assumed to carry zero flux, loaded from evotuneZeroExch.tsv. Empty cell array if that file isn't present.

See also

cmaesKcatTuning, evotuneScore

nChanged#

function n = nChanged(kcat,kcat0,varargin)

How many kcats differ from their prior by more than relTol.

The 2% default matches the threshold MATLAB's own per-source "unchanged" reporting has historically used.

Parameters

kcat, kcat0 : double: tuned and prior kcat vectors, same length.

Name-Value Arguments

relTol : double: default 0.02.

Returns

n : double: count of entries where abs(log(kcat./kcat0)) > log1p(relTol).

See also

tunePriorPenaltyWeight, parsimonyFrontier

parsimonyFrontier#

function points = parsimonyFrontier(scoreFcn,kcat,kcat0,sigma0log,varargin)

Score kcat with progressively more of it reverted to

prior.

A tuned kcat vector that fits well by rewriting every parameter is not a useful answer: changes should be few, and concentrated in the parameters least confident about. This sweeps a revert-threshold (in units of each kcat's own prior sigma) and reports how RMSE trades off against how many kcats stayed at their tuned value, so a sparser result within noise of the best fit can be preferred over the raw tuner output. See bestParsimonious to pick a point off this frontier.

Parameters

scoreFcn : function_handle: scores one kcat vector; called once per threshold. Build it against a specific ecModel/evoData, e.g. scoreFcn = @(vec) evotuneScoreVector(ecModel,tunableIdx, ... ecRxnIdsTunable,evoData,modelAdapter,params,vec); where evotuneScoreVector applies vec to ecModel.ec.kcat(tunableIdx), calls applyKcatConstraints, and returns evotuneScore's plain rmse.
kcat, kcat0, sigma0log : double: tuned vector, prior vector, and per-kcat prior sigma, all the same length.

Name-Value Arguments

thresholds : double: sigma cutoffs to sweep, ascending (default [0 1 2 3 4 6 8 10 15]). 0 reverts nothing, reproducing kcat's own score.

Returns

points : struct array: one entry per threshold, in the order given: threshold, kcat (the partially-reverted vector), nChanged, meanDev (mean movement in sigma units), rmse.

See also

bestParsimonious, nChanged

resolveTunableContext#

function [tunableIdx,ecRxnIdsTunable,kcat0,groups,sigma0log,tieMap] = resolveTunableContext(ecModel,params,tunableMask)

The tunable kcat subset plus everything derived

from it: source-group trust tiers and isozyme tie groups.

Shared by screenKcatLeverage and cmaesKcatTuning so a screen and the tuning run it feeds are always looking at exactly the same parameters.

Parameters

ecModel : struct: ecModel with a populated ec.kcat.
params : struct: modelAdapter.params.evotune (kcatSources, sigma0logSource, sigma0logDefault, tieIsozymes).
tunableMask : logical: same length as ecModel.ec.rxns; [] to consider every kcat > 0.

Returns

tunableIdx : double: indices into ecModel.ec.rxns/ec.kcat of the tunable subset.

ecRxnIdsTunable : cell :

kcat0 : double :

groups : cell: source-group name per tunable row ('unlabelled' if unmatched).
sigma0log : double: per-row prior std dev in log-space.
tieMap : double: isozymeTieMap(ecRxnIdsTunable,kcat0,'sources',groups) if params.tieIsozymes, otherwise every row is its own free parameter.

See also

screenKcatLeverage, cmaesKcatTuning, isozymeTieMap

screenKcatLeverage#

function screen = screenKcatLeverage(ecModel,varargin)

Report which kcats the data can actually speak to.

No optimisation happens here -- this only says which parameters are worth curating or tuning, and can be read on its own well before any tuning run. Each tie-group (isozyme copies sharing a prior and a source move together when params.tieIsozymes) is perturbed up and down by fold, and its leverage is the largest resulting change in RMSE. Groups are ranked by that leverage weighted by sigma0log: a trusted source needs proportionally more measured effect to rank alongside an untrusted one.

Costs one simulation per condition per tie-group probed (two probes each, up and down) plus one baseline -- comparable in scale to a full tuning run's evaluation budget, not a cheap pre-check.

Parameters

ecModel : struct: ecModel with a populated ec.kcat.

Name-Value Arguments

modelAdapter : ModelAdapter: a loaded model adapter (default: the current default model adapter).
evoData : struct: as returned by loadEvotuneData (default: loaded from modelAdapter).
tunableMask : logical: restrict which kcats are screened (default: every ec.kcat > 0).
fold : double: fold-change used to perturb each tie-group up and down (default 2.0).
nProc : double: number of parallel workers used to score candidate vectors (default 1, i.e. serial). Requires Parallel Computing Toolbox for > 1.

Returns

screen : table: one row per tie-group, sorted by rank descending: rxnId (the group's representative), nIsozymes, sourceGroup, kcat0, leverage, sigma0log, rank (leverage*sigma0log), cumLeverageShare (running total of rank as a fraction of its sum), and positions (a cell array of indices into ecModel.ec.rxns) -- what selectTunableMask needs to expand a row back into a mask.

See also

selectTunableMask, cmaesKcatTuning

selectTunableMask#

function mask = selectTunableMask(ecModel,screen,varargin)

Build a tunableMask from a screenKcatLeverage report.

Keeps the fewest highest-ranked groups whose combined leverage reaches targetImpactShare of the total -- a relative cutoff, so the same targetImpactShare selects a comparable quality of parameter set on any model, rather than an absolute leverage value or a fixed count calibrated on a different model's scale.

Pure and cheap -- no simulation, only screen is read.

Parameters

ecModel : struct :

screen : table: as returned by screenKcatLeverage.

Name-Value Arguments

targetImpactShare : double: default 0.9. Not a validated universal constant -- check the resulting mask's size against screenKcatLeverage's own curve before trusting it on a new model.

Returns

mask : logical: length numel(ecModel.ec.rxns).

See also

screenKcatLeverage, cmaesKcatTuning

sensitivityTuning#

function [model, tunedKcats] = sensitivityTuning(model, varargin)

Relax the most limiting kcats to reach a growth rate.

Iteratively increases the kcat of whichever enzyme currently limits growth the most, until the model reaches the desired growth rate, updating the tuned values in model.ec.kcat.

Parameters

model : struct: an ecModel in GECKO 3 format (with ecModel.ec structure).

Name-Value Arguments

desiredGrowthRate : double: kcats will be relaxed until this growth rate is reached (default: the experimental growth rate params.gR_exp from the model adapter).
modelAdapter : ModelAdapter: a loaded model adapter (default: the current default model adapter).
foldChange : double: kcat values will be increased by this fold-change (default 10).
protToIgnore : cell: vector of protein ids to be ignored in tuned kcats, e.g. {'P38122', 'Q99271'} (default []).
verbose : logical: whether progress should be reported (default true).

Returns

model : struct: ecModel with updated model.ec.kcat.
tunedKcats : struct: structure with information on tuned kcat values.

Notes

The tunedKcats structure has the following fields:

rxns : identifiers of reactions with tuned kcat values.
rxnNames : names of the reactions in tunedKcats.rxns.
enzymes : enzymes that catalyze the reactions in tunedKcats.rxns, whose kcat value has been tuned.
oldKcat : kcat values in the input model.
newKcat : kcat values in the output model, after tuning.
source : source label of each kcat, taken from the input model before tuning.

The model.ec.notes field will contain the original kcat value and source, unless the kcat has previously been set by sensitivityTuning, in which case the notes field remains unchanged.

Examples

% optional arguments may be given positionally or as name-value pairs:
[model, tunedKcats] = sensitivityTuning(model);
[model, tunedKcats] = sensitivityTuning(model, 'foldChange', 10);
[model, tunedKcats] = sensitivityTuning(model, desiredGrowthRate, modelAdapter, foldChange, protToIgnore, verbose);

sigmaFitter#

function [model, sigma] = sigmaFitter(model, varargin)

Fit the average enzyme saturation factor of an ecModel.

Scans candidate sigma-factors from 0.01 to 1 in steps of 0.01 and picks the one whose simulated growth rate best matches the provided experimental growth rate, then applies it to the model's protein pool constraint.

Parameters

model : struct: an ecModel in GECKO 3 format (with ecModel.ec structure).

Name-Value Arguments

growthRate : double: growth rate that should be reached (default: read from the model adapter).
Ptot : double: total cellular protein content in g/gDCW (default: read from the model adapter; if not specified there, 0.5 g/gDCW is assumed).
f : double: estimated fraction of enzymes in the model (default: read from the model adapter; if not specified there, 0.5 is assumed).
makePlot : logical: whether a plot should be made (default true).
modelAdapter : ModelAdapter: a loaded model adapter (default: the current default model adapter).

Returns

model : struct: ecModel with protein pool exchange upper bound adapted to the optimal sigma-factor.
sigma : double: optimal sigma-factor.

Examples

% optional arguments may be given positionally or as name-value pairs:
[model, sigma] = sigmaFitter(model);
[model, sigma] = sigmaFitter(model, 'makePlot', false);
[model, sigma] = sigmaFitter(model, growthRate, Ptot, f, makePlot, modelAdapter);

truncateValues#

function table = truncateValues(table,cols)

Round values in selected columns to six significant digits.

Rounds each value in the given columns of a cell array or table to six significant digits, based on the value's order of magnitude.

Parameters

table : cell or table: a cell array or table where some columns have values that should be truncated.
cols : double: index or indices of columns with values to be truncated.

Returns

table : cell or table: the input with the selected columns rounded.

tunePriorPenaltyWeight#

function sweepTable = tunePriorPenaltyWeight(ecModel,varargin)

Sweep priorPenaltyWeight and report the

fit/reproducibility trade-off, so the value doesn't have to be picked blind.

priorPenaltyWeight charges for moving a kcat away from its prior, and mainly exists to keep large corrections reproducible across seeds rather than landing on an arbitrary point along a flat direction. How strong a charge is enough is model- and data-specific: too little and the search's largest corrections disagree across seeds; too much and the search stops finding real corrections at all. This runs cmaesKcatTuning at every value in candidates, each at every seed in seeds, all against the same tunable set, and reports both fit and cross-seed reproducibility per value.

This is exactly as expensive as it sounds: numel(candidates)*numel(seeds) full cmaesKcatTuning runs, none of them reusable across candidates. Pass a screen you already have to at least avoid recomputing that part.

Parameters

ecModel : struct: ecModel with a populated ec.kcat. Not mutated (MATLAB structs are copied on assignment; every candidate/seed pair starts fresh from this input regardless of what earlier runs did).

Name-Value Arguments

modelAdapter, evoData, tunableMask, screen, targetImpactShare, popsize, :

nProc, verbose: as in cmaesKcatTuning.
candidates : double: priorPenaltyWeight values to try (default [0.0, 0.01, 0.03, 0.1]).
seeds : double: seeds to run each candidate at (default [0, 1]). Reproducibility columns pool every pair of seeds, so more than two sharpens them but costs proportionally more.
foldThreshold : double: a kcat counts as "moved" past this fold change from its prior (default 2.0). Matches nChanged's convention loosely (a different quantity, same spirit).

Returns

sweepTable : table: one row per candidate: priorPenaltyWeight, nSeeds, distanceMean, distanceSpread (max minus min across seeds), nChangedMean, nMovers (kcats past foldThreshold in either seed of a pair, summed over every pair), nBothMoved (past it in both), pctDirectionAgree, medianFoldSpread and maxFoldSpread (among movers, the largest fold-change between what two seeds landed on), and fitCost (distanceMean relative to the best candidate's, as a fraction).

See also

cmaesKcatTuning, screenKcatLeverage, selectTunableMask, nChanged