Tuning against experimental data#

Note

Python code, with a MATLAB equivalent Every code block on this page is Python. GECKO MATLAB has an equivalent module, src/kcat_tuning/evotune/, with the corresponding functions loadEvotuneData, screenKcatLeverage, selectTunableMask, cmaesKcatTuning, tunePriorPenaltyWeight and, for the assignment-review step below, reviewKcatAssignment; this page does not cover the MATLAB call signatures.

Growth-rate tuning raises kcats one at a time until a single target growth rate is reached. cmaes_kcat_tuning instead fits many kcats at once, with CMA-ES, against a full set of measured growth rates and fluxes across several conditions, weighted by how much each kcat's source is trusted. This page covers one recommended path through that fitting process end to end, in geckopy. It is not a full reference; the docstrings in geckopy.kcat_tuning.evotune document every option.

Functions on this page#

Python
`screen_kcat_leverage`	report which kcats the data can inform, before any tuning runs
`review_assignment`	check kcat values for assignment mistakes, independent of the data
`cmaes_kcat_tuning`	search the tunable kcats with CMA-ES for the best fit to measured data
`result.corrections`	turn the tuning result into a reviewable list of changes

What this does#

An ecModel's kcats mostly come from databases (BRENDA), predictors (DLKcat, OpenKineticsPredictor), or manual curation (see Gathering kcats and Applying kcats), and they do not on their own reproduce measured growth. Three functions, run in sequence, turn measured data into corrections:

screen_kcat_leverage reports which kcats the data can actually inform: no optimization runs yet, only "if this kcat were different, how much would the fit change." This is useful for curation on its own, before any tuning run.
select_tunable_mask turns that report into the set of kcats a tuning run searches over. It runs automatically if skipped.
cmaes_kcat_tuning searches that set with CMA-ES for the kcat vector that best matches the data, weighted by how much each kcat's source is trusted.

A fourth function, review_assignment, is complementary rather than part of that sequence: it checks the kcats themselves for assignment mistakes, independent of any measured data, since the tuner's largest corrections are usually places the assignment went wrong rather than discoveries about biology (see Check the assignments behind the top kcats below).

The result is usually a short list of corrections, each with a before/after value and a source, not a wholesale rewrite of the model. A run that changes nearly every kcat by a little has not found anything; a run that changes a handful by a lot, concentrated in the least-trusted sources, has.

Before you start#

A working ecModel and ModelAdapter project: an existing model_adapter.toml and models/*.yml (see Getting started and Building an empty ecModel).
pip install geckopy[evotune]. The CMA-ES search (cma) is an optional dependency, not part of the base install.
Experimental data: at least one of a set of measured exchange fluxes per condition, or a set of measured maximum growth rates per carbon source. Both together works and is preferred, since flux data pins the shape of metabolism and growth rate pins its rate.

Prepare the experimental data#

Three tab-separated files go under <project>/data/. All three are optional; a missing file is skipped rather than an error, but at least one of the first two is required.

File	Contents	Required
`evotuneFluxData.tsv`	Measured exchange fluxes, one row per condition.	At least one of these two
`evotuneMaxGrowth.tsv`	Measured maximum growth rate, one row per condition, one active carbon source per row (at `-1000`, unconstrained).
`evotuneZeroExch.tsv`	Reaction IDs assumed to carry zero flux in every condition.	No

The two flux files share one column layout, the same parser geckopy uses for regular flux data:

Condition   Ptot   grRate   glucose (r_1714)   fructose (r_1709)   ...   evotuneRMSEweight   source
glucose     NaN    0.41     -1000              NaN                ...   1                     DLKcat
fructose    NaN    0.338    NaN                -1000              ...   1                     DLKcat

Condition names the row; grRate is the measured growth rate; Ptot is measured protein content (NaN if not measured); each <met> (<rxn>) column is a measured or fixed flux for that exchange reaction; evotuneRMSEweight scales how much that row counts toward the RMSE (usually 1); and source is free text carried through for bookkeeping. evotuneZeroExch.tsv is simpler: a Rxns header, then one reaction ID per line.

Configure evotune in model_adapter.toml#

Every field has a default, so this section can be omitted entirely to start. A configured example, in plain terms:

[evotune]
sigma0_log_default = 0.3   # trust for any source not listed below
max_growth_weight = 2.0    # weight the growth-rate data double against flux data

[evotune.source_groups.okp]
sources = ["OpenKineticsPredictor"]
match_okp = true

[evotune.source_groups.brenda]
sources = ["brenda"]

[evotune.source_groups.custom]
sources = ["custom"]

[evotune.sigma0_log_source]
okp = 0.25
brenda = 0.2
custom = 0.1

What each setting controls:

Trust tiers (source_groups, sigma0_log_source and sigma0_log_default). Group kcats by how much their origin is trusted, and give each group a standard deviation in log-space: a smaller value means more trust, so the search needs stronger evidence before moving that kcat far. custom (curated by hand) at 0.1 is three times more trusted than an unlabelled prediction at 0.3. This is the one section with no universal default; the right tiers depend on the model's own source labels.
max_growth_weight (default 1.0). The combined score is (rmse_flux + w * rmse_max_growth) / (w + 1). At 2, growth-rate error counts double against flux error, useful when growth rate is measured across many carbon sources but flux is mostly measured for one.
prior_penalty_weight (default 0.03, already on). Adds w * mean(((log k - log k0) / sigma0)^2) to the score: a charge for moving a kcat, scaled by how much its prior is trusted. Without this term the fit is often flat along many directions: several very different kcat vectors fit equally well, and the search returns an arbitrary point along that flat direction, one that does not reproduce when rerun with a different seed. The default is calibrated so that large corrections agree in direction and land within a handful of fold of each other across seeds, at a small cost in raw fit. Setting it to 0 scores on fit alone, which is also what a run reproducing GECKO MATLAB, which has no such term, must do. The right strength depends on the model and the data, since it scales with sigma0_log and the objective itself; tune_prior_penalty_weight runs a candidate sweep and reports fit cost against cross-seed reproducibility so the value does not have to be picked blind (see its docstring). It costs one full tuning run per candidate per seed, so treat it as an occasional calibration step rather than part of routine tuning.
tie_isozymes (default true). Several ecModel reactions are one enzymatic reaction split across isozyme copies. When copies share a prior value and a source, no experimental condition can distinguish between them, so an untied search can assign them a difference that carries no biological meaning. Tying gives such copies one shared kcat.
max_generations / rmse_threshold. This search's stopping conditions: run at most this many CMA-ES generations, or stop earlier once the best RMSE reaches rmse_threshold (negative never stops early).

Screen which kcats the data can inform#

from geckopy import ModelAdapter, load_ec_model
from geckopy.kcat_tuning.evotune import screen_kcat_leverage

adapter = ModelAdapter.from_folder("path/to/project")
model = load_ec_model(adapter=adapter)  # models/ecModel.yml by default

screen = screen_kcat_leverage(model, adapter=adapter, n_proc=8)
screen.drop(columns="_positions").head(20)

This perturbs every tunable kcat, moving isozyme copies together per tie_isozymes, up and down, and measures how much the fit changes; no tuning happens yet. The result is a table, one row per kcat or tied group, ranked by that leverage weighted by trust: the reaction ID, how many isozyme copies it covers, its source group, its current value, its leverage, and cum_leverage_share, the running fraction of total leverage carried by this row and every row above it. The top rows are what a curator should look at first, independent of whether a tuning search ever runs: on ecYeastGEM, a handful of kcats routinely carry the majority of what any tuning run could achieve, and they tend to be implausible database values rather than genuine biology.

This costs one simulation per condition per kcat group, roughly comparable in scale to a full tuning run's budget; expect tens of minutes on a genome-scale model, not a quick check. n_proc parallelizes it the same way tuning does.

Check the assignments behind the top kcats#

Before spending a tuning run on a kcat, check whether the value was assigned correctly in the first place. review_assignment flags source-level problems: a value with no kcat row behind it for its EC, one number reused across many unrelated reactions, an EC's maximum standing in for a whole distribution, an implausible magnitude, or a curated value that repeats what the database already said. These are corrections to the source, so they belong in whatever database or fuzzy-matching step in Gathering kcats produced the value; fixing them there propagates to every model built from the same data, unlike a value the tuner moves for this one model alone.

from geckopy.gather_kcats.review import EcStats, review_assignment, findings_tsv

# One EcStats per EC code queried while gathering kcats, summarising
# whatever kcat database (usually BRENDA) the assignment came from.
ec_stats = {
    "1.1.1.1": EcStats(n_kcat=12, kcat_max=340.0, kcat_median=45.0,
                       values=frozenset({340.0, 45.0, 12.0}), n_sa=0),
    # ...
}

findings = review_assignment(
    model.ec.rxns, model.ec.kcat, model.ec.source, model.ec.eccodes,
    ec_stats=ec_stats,
)
print(findings_tsv(findings))

Each Finding names a reaction, its checks (no-ec-evidence, repeated-value, ec-maximum, magnitude, custom-duplicate), and a human-readable detail. Isozyme copies the assignment could not distinguish, the same grouping tie_isozymes uses above, are reported once, not once per copy. Standard fallbacks and tied groups are deliberate choices rather than mistakes, so neither is reported unless explicitly requested via include.

Findings rank by leverage when it is passed in, reusing the screen table from the previous step, expanded from tie-groups back to individual model.ec.rxns positions via its _positions column, so that looking odd and actually mattering are not confused; without it, rows come back in input order, and looking odd is all that remains. coverage (a share of the flagged rows' total leverage) or top (a row count) then truncates the report to a length worth a curator's time: coverage adapts to the model, keeping however many rows it takes to carry that share of the flagged leverage, where a fixed top does not. On ecYeastGEM the underlying checks flag around 1,900 of 4,834 kcats before any truncation, which on its own is still too long a report to read; tightening the checks themselves does not fix this nearly as well as truncating by leverage does, since looking odd (what the checks measure) and actually mattering (what leverage measures) are unrelated properties, so a stricter check threshold drops high-leverage entries as readily as low-leverage ones.

Tune#

from geckopy import save_ec_model
from geckopy.kcat_tuning.evotune import cmaes_kcat_tuning

result = cmaes_kcat_tuning(
    model, adapter=adapter, screen=screen, n_proc=8, seed=0,
)

save_ec_model(model, "ecModel_tuned.yml", adapter=adapter)

Passing the screen from the earlier step avoids recomputing it; omitting screen entirely makes cmaes_kcat_tuning compute one itself. Either way, the tunable set is built by select_tunable_mask, which keeps the fewest highest-ranked kcats whose combined leverage reaches target_impact_share (default 0.9) of the total: a cutoff relative to the model's own achievable improvement, not a fixed count or an absolute leverage value, so the same setting means a comparable thing on a different model. Pass tunable_mask directly to fix the set by hand instead.

model is mutated in place: on return, every tunable kcat carries the best value CMA-ES found, and the model's kcat constraints are already applied, so nothing further is needed before simulating or saving it. Beyond target_impact_share, the only other knobs are popsize (CMA-ES's population size, defaulting to its own dimension-scaled choice rather than a value tuned for one particular model), n_proc, and seed; everything else comes from the same configuration as the previous step, so there is nothing new to configure once that section is read.

If growth needs special handling for the organism, for example forcing anaerobic conditions or scaling the protein pool with a condition-specific biomass composition, pass make_anaerobic and/or change_protein_biomass callables; geckopy has no organism-agnostic default for either. tutorials/full_ecModel/code/anaerobic.py in the geckopy repository has a worked example.

Read the result#

result is an EvotuneResult:

rmse_trace / objective_trace: best-so-far plain-fit RMSE and best-so-far optimized objective, per generation (identical unless prior_penalty_weight is nonzero). A flat trace for many generations before the run ends means it stalled, not converged.
new_kcat / old_kcat / rxns / groups: the tuned and prior kcats over the searched set, parallel arrays, with each kcat's trust tier. Tied isozyme copies share one value in new_kcat.
converged: whether rmse_threshold was reached, or max_generations was hit first.
n_generations: how many generations actually ran.

Before trusting a tuned kcat, check more than the final RMSE:

How many kcats changed, and by how much. A result that moves nearly everything by a little has not identified anything. geckopy.kcat_tuning.evotune.parsimony has the tools for this: n_changed, fold_change, source_movement.
Impact share. parsimony.impact_share reports what fraction of the total achievable improvement the changed kcats actually carry, distinguishing "few, large and consequential" from mere sparsity. This is a different, unweighted quantity from screen_kcat_leverage's trust-weighted cum_leverage_share.
Reproducibility. Rerun with a different seed. Large corrections (beyond roughly two-fold) should land in the same direction and within a small factor of each other; if they do not, the search has found a flat direction rather than a real correction, and prior_penalty_weight is the first setting to reach for.
Whether it generalizes. If data exists that was not fit on, score the tuned model against it and compare to the untuned model. A large improvement on fitted conditions with no improvement elsewhere is a sign of overfitting to a thin dataset, not a fixed model.

new_kcat/old_kcat is thousands of numbers, most of them unchanged, not something a reviewer can read directly. result.corrections() turns it into the actual review artifact: one Correction row per changed kcat (reaction, EC code, source, prior/tuned value, fold change), ranked by leverage rather than by how far a kcat moved, so parameters the data cannot see, which are free to drift furthest and are the least reliable, sort to the bottom instead of the top:

from geckopy.kcat_tuning import annotate_from_model, corrections_tsv

# screen is the table from "Screen which kcats the data can inform" above; reuse
# it rather than recomputing leverage.
leverage_by_rxn = dict(zip(screen["rxn_id"], screen["leverage"]))
leverage = [leverage_by_rxn.get(r, 0.0) for r in result.rxns]

annotations = annotate_from_model(model, result.rxns)
rows = result.corrections(
    leverage=leverage, names=annotations["names"], ec_codes=annotations["ec_codes"],
)
open("corrections.tsv", "w").write(corrections_tsv(rows))

cumulative_share on each row is the running fraction of total leverage the list has accounted for so far, so a reader can stop reading once it stops climbing. This is the tuning-time counterpart to review_assignment above: that step flags kcats that need curation before a run; this one reports which ones the run actually changed, and how much each change can be trusted.

Why CMA-ES, not ABC-SMC#

DLKcat's own Bayesian approach to kcat tuning implements ABC-SMC (Approximate Bayesian Computation, Sequential Monte Carlo); see SysBioChalmers/DLKcat/BayesianApproach. geckopy's kcat_tuning.evotune module ported that approach faithfully and tested it side by side with CMA-ES, a plain optimizer, on the same objective, the same screened parameter set and the same experimental data (ecYeastGEM, 41 conditions). CMA-ES won on every criterion checked, which is why it is what geckopy ships.

ABC-SMC is a sampler: each generation it draws a batch of candidate kcat vectors from a proposal distribution, scores them, and keeps the best fraction to seed the next generation's proposals. That machinery, a posterior approximation, a sample schedule and truncation selection, earns its keep when the goal is a distribution over plausible parameter values and the model is cheap enough to sample densely across the whole parameter space. Once a screen has reduced the problem to the kcats the data can actually inform (a few dozen to a few hundred, out of thousands), that premise stops holding: there is no longer a broad space to sample across, and each evaluation is expensive, one FBA solve per condition. At that point the problem is a direct search for the best-fitting vector, which is what an optimizer is for and not what ABC-SMC is for.

On the same model, the same 112-parameter screen, the same objective and three seeds each:

	Distance	Spread across seeds
ABC-SMC	0.9156	+/- 0.0549
CMA-ES	0.7974	+/- 0.0048

CMA-ES fits better by 3.7 standard errors, with a spread more than ten times tighter. The gap is not only fit quality; it also determines whether an individual tuned kcat means anything. Both methods leave most kcats on flat directions the data cannot pin down, but ABC-SMC's posterior sampling spreads its samples across the whole flat region: across three ABC-SMC seeds, lanosterol synthase, the single highest-leverage kcat in the model, took values of 0.338, 20.3 and 0.0781 s^-1, a 260-fold spread, at distances differing by only 0.009. A number that unstable cannot be reported as a finding.

CMA-ES does not remove flat directions; nothing can, since they are a property of the data rather than of the search method. A direct optimizer converges to one point per flat direction instead of sampling across it, and prior_penalty_weight, a Tikhonov-style penalty on moving away from the prior, weighted by how much each source is trusted, picks out which point.