kcat sensitivity tuning (Python)#

geckopy objects in geckopy.kcat_tuning.

Functions#

Function	Summary
`annotate_from_model`	Names, EC codes and sources for `rxn_ids`, read off an ecModel.
`apply_ties`	Give every tied position its representative's value.
`base_reaction`	The reaction an `ec.rxns` entry belongs to, suffix removed.
`best_parsimonious`	The fewest-changes point whose RMSE is within `tolerance` of the
`build_sigma0_log`	Per-row prior std dev in log-space, from source-group membership.
`classify_kcat_source`	Map one `ec.source` string to a trust-tier group name.
`classify_kcat_sources`	Vectorised `classify_kcat_source` over `ec.source`.
`cmaes_kcat_tuning`	Tune kcats against experimental data with CMA-ES.
`compute_excarbon`	Per-reaction carbon-weight lookup for RMSE weighting.
`ConditionSimResult`	One condition's simulated outcome.
`Correction`	One kcat the tuning changed, with its provenance.
`corrections`	The changed kcats, most consequential first.
`corrections_tsv`	The corrections as a TSV document, header included.
`dataset_rmse`	RMSE (mean, per-condition array) for one tuning dataset.
`evotune_distance`	Combine both datasets' RMSE into one score, per `abc_max.m`.
`EvotuneData`	Experimental data used by evotune kcat tuning.
`EvotuneResult`	Outcome of a `cmaes_kcat_tuning` run.
`find_max_value`	Find the maximum kinetic value across one or more EC numbers.
`fit_sigma`	Sweep sigma and pick the one matching `growth_rate` best.
`fold_change`	Per-kcat fold change from the prior, always >= 1.
`FrontierPoint`	One (parsimony, fit) trade-off, from `parsimony_frontier`.
`identifiability_mask`	Which kcats have earned the right to be tuned.
`impact_share`	Fraction of the total screened impact carried by changed kcats.
`isozyme_tie_map`	Index of the parameter each position follows.
`kcat_bounds`	Biologically plausible bounds for proposed kcats, in 1/s.
`load_evotune_data`	Load the three evotune-kcat-tuning experimental data files.
`MaskSizePoint`	Tuned fit at one candidate `target_impact_share`.
`movement_in_sigma`	Per-kcat displacement from the prior, in prior standard deviations.
`n_changed`	How many kcats differ from their prior by more than `rel_tol`.
`n_free`	How many parameters remain free once ties are applied.
`parsimony_frontier`	Score `kcat` with progressively more of it reverted to prior.
`recommend_target_impact_share`	Pick the smallest share past which tuning stops paying for itself.
`revert_below`	Put every kcat that moved less than `threshold` sigma back on
`screen_kcat_leverage`	Report which kcats the data can actually speak to.
`select_tunable_mask`	Build a `tunable_mask` from a `screen_kcat_leverage` report.
`sensitivity_tuning`	Iteratively bump the most-limiting kcat by `fold_change` until
`SigmaFitterResult`	Outcome of a sigma scan.
`simulate_evotune_dataset`	Simulate every condition in one tuning dataset.
`source_movement`	Median fold change and near-prior fraction, per source group.
`sweep_tunable_mask_size`	Tune at each of `target_impact_shares` and report the fit.
`truncate_values`	Round values in `arr` to ~6 significant figures.
`tune_prior_penalty_weight`	Sweep `prior_penalty_weight` and report the fit/reproducibility
`TunedKcatsResult`	Outcome of a sensitivity-tuning run.

Reference#

annotate_from_model#

annotate_from_model(model, rxn_ids: Sequence[str])

Names, EC codes and sources for rxn_ids, read off an ecModel.

An isozyme copy carries a _EXP_n suffix that the underlying reaction does not, so the name is looked up under the stripped id. Anything absent is reported as "".

apply_ties#

apply_ties(particles: np.ndarray, tie_map: np.ndarray)

Give every tied position its representative's value.

particles is (n_params,) or (n_params, n_particles); the array is modified in place and returned, so this can sit directly in a sampling loop.

base_reaction#

base_reaction(rxn_id: str)

The reaction an ec.rxns entry belongs to, suffix removed.

best_parsimonious#

best_parsimonious(points: Sequence[FrontierPoint], tolerance: float=0.02)

The fewest-changes point whose RMSE is within tolerance of the

best on the frontier.

Defaults to 2%, which is the scale of the alternate-optimum noise measured on this model -- below that, two vectors are not meaningfully different in fit, so the sparser one wins.

build_sigma0_log#

build_sigma0_log(groups: np.ndarray, params: EvotuneParams)

Per-row prior std dev in log-space, from source-group membership.

Mirrors MATLAB's sigma0log = sigma0logDefault * ones(...); sigma0log(uniqKcatParams) = sigma0logSource(kcatSourceIdx(...)): every row starts at params.sigma0_log_default, then rows whose group is a real (non-UNLABELLED_GROUP) entry get overridden by that group's sigma0_log_source value.

Parameters

groups: Per-row group names from classify_kcat_sources.
params: Supplies sigma0_log_default and sigma0_log_source.

Returns

numpy.ndarray of float, same shape as ``groups``. :

classify_kcat_source#

classify_kcat_source(source: str, params: EvotuneParams, *, okp_method: Optional[str]=None)

Map one ec.source string to a trust-tier group name.

Matching is case-insensitive, mirroring MATLAB's strcmpi(ecModel.ec.source, kcatSources{i}), and ignores the bracketed fuzzy-match detail geckopy appends to a database source. ec.source reads "brenda (wc=0, origin=1)" where MATLAB's reads "brenda", so matching the whole string would leave every fuzzy-matched kcat unlabelled and silently on the default prior width -- which is most of a model's kcats. A source_groups entry listing the full string still matches it.

Parameters

source: One ec.source entry.
params: Supplies source_groups.
okp_method: The project's configured OpenKineticsPredictor method (OkpParams.method), if any. A group with match_okp=True also matches source when it equals this (case-insensitive).

Returns

str: A key of params.source_groups, or UNLABELLED_GROUP if nothing matches -- callers then fall back to params.sigma0_log_default / shrink_thr_default / force_prior_thr_default.

classify_kcat_sources#

classify_kcat_sources(ec_sources: list[str], params: EvotuneParams, *, okp_method: Optional[str]=None)

Vectorised classify_kcat_source over ec.source.

Returns

`numpy.ndarray of str, shape ``(len(ec_sources),)``` :

cmaes_kcat_tuning#

cmaes_kcat_tuning(model: 'EcModel', *, adapter: Optional['ModelAdapter']=None, params: Optional['EvotuneParams']=None, evotune_data: Optional[EvotuneData]=None, okp_method: Optional[str]=None, bio_rxn: Optional[str]=None, make_anaerobic: Optional[Callable[['EcModel'], None]]=None, change_protein_biomass: Optional[Callable[['EcModel', float], None]]=None, tunable_mask: Optional[np.ndarray]=None, screen: Optional[pd.DataFrame]=None, target_impact_share: float=0.9, popsize: Optional[int]=None, n_proc: Optional[int]=None, seed: Optional[int]=None, verbose: bool=True)

Tune kcats against experimental data with CMA-ES.

Configuration comes entirely from EvotuneParams: trust tiers, max_growth_weight, prior_penalty_weight and tie_isozymes control the search, and max_generations/rmse_threshold are its stopping conditions.

Parameters

model: EcModel with a populated ec.kcat. Mutated in place: on return, tunable rows carry the best kcat vector CMA-ES found, with kcat constraints already applied.
tunable_mask, screen, target_impact_share: Control which kcats are searched over, in order of precedence. Pass tunable_mask to fix the set yourself. Otherwise, a mask is built by select_tunable_mask at target_impact_share, from screen if given (so a previously computed screen_kcat_leverage report need not be recomputed) or from a fresh screen otherwise.
popsize: CMA-ES population size. Defaults to cma's own dimension-scaled default (4 + floor(3 ln n)), so it adapts to however many free parameters the mask/tying produced rather than a value calibrated on one particular model.
n_proc: Number of worker processes for scoring each generation's candidates. Defaults to cobra.Configuration().processes. 1 runs the original serial path (no Pool at all). See the module docstring's "Parallel scoring" section.
seed: If given, seeds the search, for reproducible runs.
verbose: Whether per-generation progress is logged at INFO.
make_anaerobic, change_protein_biomass: Forwarded to simulate.simulate_evotune_dataset -- see its docstring; geckopy has no generic organism-agnostic implementation of these yet. See the module docstring's "Parallel scoring" section for a picklability caveat when n_proc>1.

Returns

EvotuneResult: rxns/old_kcat/new_kcat/groups cover the selected tunable set, tied groups included (their members share one value in new_kcat). rmse_trace/objective_trace record the best-so-far plain RMSE and optimised objective per generation.

Raises

ValueError: 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.

compute_excarbon#

compute_excarbon(model: 'EcModel', rxn_ids: Iterable[str], bio_rxn_id: str)

Per-reaction carbon-weight lookup for RMSE weighting.

The biomass reaction gets BIOMASS_CARBON_EQUIV. Every other reaction gets the carbon count of its single exchanged metabolite (via cobrapy's parsed Metabolite.elements, not a hand-rolled formula parser). A missing/unparseable formula, or a computed carbon count of exactly 0 (e.g. O2, water), both fall back to 1 -- matching MATLAB's Ematrix(isnan(Ematrix)) = 1 and excarbon(excarbon == 0) = 1 respectively.

Parameters

model: Model providing the metabolite formulas.
rxn_ids: Reaction IDs to compute weights for (typically every reaction referenced by a tuning dataset: its exchange columns, its zero-flux list, and the biomass reaction).
bio_rxn_id: The biomass reaction ID.

Returns

dict[str, float] :

ConditionSimResult#

One condition's simulated outcome.

feasible=False (LP infeasible/unbounded, or no biomass flux) means growth/exch_fluxes/block_fluxes are meaningless; .distance applies the fixed infeasibility penalty for that row instead of reading them.

Correction#

One kcat the tuning changed, with its provenance.

corrections#

corrections(kcat: np.ndarray, kcat0: np.ndarray, rxn_ids: Sequence[str], *, sources: Optional[Sequence[str]]=None, names: Optional[Sequence[str]]=None, ec_codes: Optional[Sequence[str]]=None, leverage: Optional[np.ndarray]=None, rel_tol: float=0.02)

The changed kcats, most consequential first.

Parameters

kcat, kcat0: Tuned and prior kcat vectors, same length as rxn_ids.
sources, names, ec_codes: Per-kcat provenance. Missing entries are reported as "" rather than dropped, so a row never silently disappears.
leverage: Per-kcat effect on the distance from a one-at-a-time screen (see .tuning.screen_kcat_leverage). If given, rows are ordered by it and the share columns are filled; otherwise rows are ordered by fold change and the shares are 0.
rel_tol: A kcat within 1 + rel_tol fold of its prior is unchanged.

Returns

list of :class:Correction, ordered most consequential first. :

corrections_tsv#

corrections_tsv(rows: Sequence[Correction])

The corrections as a TSV document, header included.

dataset_rmse#

dataset_rmse(flux_data: 'FluxData', sims: list[ConditionSimResult], *, constrain: bool, excarbon: dict[str, float], bio_rxn_id: str, penalty: float=INFEASIBLE_PENALTY)

RMSE (mean, per-condition array) for one tuning dataset.

Parameters

flux_data: The dataset being scored (evotune_data.flux_data or evotune_data.max_grate).
sims: One .simulate.ConditionSimResult per flux_data.conds row, e.g. from .simulate.simulate_evotune_dataset.
constrain: Must match the value passed to .simulate.simulate_evotune_dataset for these sims: True includes carbon-weighted exchange-flux terms (flux data), False scores growth alone (max-growth data).
excarbon: Carbon-weight lookup from compute_excarbon, covering every reaction referenced by sims' exch_fluxes / block_fluxes plus bio_rxn_id. A reaction missing from the lookup falls back to a weight of 1.
bio_rxn_id: The biomass reaction ID (for reading flux_data.gr_rate's carbon-equivalent counterpart).
penalty: RMSE substituted for an infeasible condition.

Returns

rmse : float: Weighted (by flux_data.evotune_rmse_weight, if present) mean of the per-condition RMSEs.
rmse_list : numpy.ndarray: The per-condition RMSEs (post-penalty, pre-weighting), shape (len(flux_data.conds),).

evotune_distance#

evotune_distance(evotune_data: 'EvotuneData', *, flux_sims: list[ConditionSimResult] | None, max_grate_sims: list[ConditionSimResult] | None, excarbon: dict[str, float], bio_rxn_id: str, penalty: float=INFEASIBLE_PENALTY, max_growth_weight: float=1.0)

Combine both datasets' RMSE into one score, per abc_max.m.

flux_sims must be given (non-None) iff evotune_data.flux_data is not None, and likewise for max_grate_sims / evotune_data.max_grate. A dataset absent from evotune_data doesn't contribute to the combined RMSE (MATLAB: a missing dataset's rmse_* stays [] and drops out of validIdx).

Parameters

max_growth_weight

Relative weight of the max-growth dataset against the flux dataset, giving (rmse_flux + w * rmse_max_grate) / (w + 1). At 2 the eight max-growth conditions carry twice the weight of the 33 flux conditions, which is the intended emphasis here: the max-growth set is where the carbon-source diversity lives, while 30 of the 33 flux conditions are glucose.

This is deliberately not MATLAB's convention. abc_max.m applies weights = [maxGrowthWeight, 1] to [rmse_flux, rmse_maxGrate], so there the same field scales the flux term. A run reproducing MATLAB must pass 0.5 here to match its maxGrowthWeight = 2.

Returns

rmse : float: Weighted mean of whichever dataset RMSEs are present. NaN if neither dataset is present.
detail : dict[str, numpy.ndarray]: Per-condition RMSE arrays, keyed "flux_data" / "max_grate" for whichever datasets were scored.

EvotuneData#

Experimental data used by evotune kcat tuning.

Attributes

flux_data: Per-condition growth/exchange-flux measurements, from evotuneFluxData.tsv. None if that file isn't present.
max_grate: Per-condition maximum-growth-rate measurements (one active carbon source per condition, at -1000 i.e. unconstrained), from evotuneMaxGrowth.tsv. None if that file isn't present.
zero_flux: Reaction IDs assumed to carry zero flux, from evotuneZeroExch.tsv. Empty if that file isn't present.

EvotuneResult#

Outcome of a cmaes_kcat_tuning run.

The model is mutated in place (the best kcat vector found is applied to model.ec.kcat), and this result records the run's history for inspection.

Attributes

rxns: Tunable ec.rxns IDs, in the column order used throughout (rxns[i] corresponds to old_kcat[i]/new_kcat[i]/ groups[i]). Tied isozyme copies are all listed, sharing one value in new_kcat.
old_kcat, new_kcat: kcat values before/after tuning, parallel to rxns.
groups: Source-group name per tunable row (from priors.classify_kcat_sources).
rmse_trace: Best-so-far plain RMSE, per generation. Always the plain RMSE, so runs stay comparable across penalties and against MATLAB even when the search used a penalised objective.
objective_trace: Best-so-far value of the quantity the search actually minimised, per generation. Identical to rmse_trace when params.prior_penalty_weight is 0.
n_generations: Number of generations actually run.
converged: True if rmse_trace[-1] <= params.rmse_threshold; False if max_generations was reached first.

find_max_value#

find_max_value(ec_string: str, brenda: 'BrendaData')

Find the maximum kinetic value across one or more EC numbers.

Searches both brenda.kcat_max and brenda.sa_max for each EC token in ec_string and returns the overall maximum together with its organism and parameter type. Wildcards in EC tokens (-) trigger prefix matching: "1.1.1.-" matches any code starting with "1.1.1.".

Ported from GECKO MATLAB: src/geckomat/kcat_sensitivity_analysis/findMaxValue.m.

Parameters

ec_string: Space-separated EC tokens. Each token may optionally start with "EC" (matching the MATLAB convention) and may contain - wildcards.
brenda: Pre-loaded BRENDA tables.

Returns

value: The overall maximum kinetic value (1/s).
organism: The organism contributing the maximum.
parameter: Either "K_cat" (came from the kcat table) or "SA*Mw" (came from the SA-derived kcat table).

If no EC matches anything in either table, returns :

```(0.0, "", "")``.` :

fit_sigma#

fit_sigma(model: 'EcModel', *, growth_rate: Optional[float]=None, p_tot: Optional[float]=None, f: Optional[float]=None, n_sigma_steps: int=100, method: Literal['grid', 'bisect']='bisect')

Sweep sigma and pick the one matching growth_rate best.

Growth is monotone non-decreasing in sigma, so method="bisect" (default) finds the best sigma in about log2(n) solves; the returned grids then contain only the sigmas actually evaluated. method="grid" keeps the legacy [1/n, 2/n, ..., 1.0] sweep for diagnostic plotting (every grid point evaluated).

The best sigma is re-applied to model before return.

Ported from GECKO MATLAB: src/geckomat/kcat_sensitivity_analysis/sigmaFitter.m.

Parameters

model: EcModel with the protein pool machinery installed and the objective set (typically the biomass reaction). Mutated in place: the optimal sigma is applied via set_prot_pool_size before return.
growth_rate: Experimental growth rate to match. Defaults to model.adapter.params.gr_exp.
p_tot: Total protein content (g/gDCW). Defaults to model.adapter.params.p_tot.
f: Mass fraction of model enzymes. Defaults to model.adapter.params.f.
n_sigma_steps: For method="grid", the number of sigma values tried in (0, 1] (i / n for i = 1..n). For method="bisect", it caps the number of bisection iterations (the default 100 is far more than the ~20 needed to converge). Default 100.
method: "bisect" (default) finds the best sigma in ~log2(n) solves using the monotonicity of growth in sigma; "grid" sweeps the full grid.

Returns

SigmaFitterResult :

Raises

ValueError: If model.adapter is None and any default is needed; if n_sigma_steps is < 1; or if method is unknown.

fold_change#

fold_change(kcat: np.ndarray, kcat0: np.ndarray)

Per-kcat fold change from the prior, always >= 1.

max(k/k0, k0/k). Reported in preference to movement_in_sigma because sigma0_log differs by source (0.10 for custom up to 0.30 for unlabelled), so dividing by it hides the per-source behaviour: a vector that moves every source by the same number of sigmas has moved the trusted ones three times less in fold terms. Only the ranking of the sigmas across sources is meaningful, not their absolute values.

FrontierPoint#

One (parsimony, fit) trade-off, from parsimony_frontier.

identifiability_mask#

identifiability_mask(drmse: np.ndarray, sigma0_log: np.ndarray, threshold: float)

Which kcats have earned the right to be tuned.

Keeps a kcat when |dRMSE| * sigma0_log > threshold, where dRMSE is its one-at-a-time effect on the distance. Multiplying by sigma0_log makes a trusted source clear a proportionally higher bar before it may move: with custom at 0.10 and unlabelled at 0.30, a custom kcat needs three times the measured effect of an unlabelled one to qualify.

This uses only the ranking of the sigmas across sources, which is the part that carries meaning -- their absolute values are a convention. One global threshold sets the overall strictness; the per-source differentiation is derived, not tuned.

isozyme_tie_map#

isozyme_tie_map(rxn_ids: Sequence[str], kcat0: np.ndarray, sources: Optional[Sequence[str]]=None, *, rel_tol: float=1e-09)

Index of the parameter each position follows.

Positions that are their own representative are free; the rest follow one. Copies group together when they share a base reaction, a prior value (to rel_tol) and -- when sources is given -- a source group.

Returns

`numpy.ndarray of int, one entry per position. ``tie_map[i] == i``` :

for a free parameter, otherwise the index it follows. The :

representative of a group is its lowest index, so the map is :

idempotent: ``tie_map[tie_map] == tie_map``. :

kcat_bounds#

kcat_bounds(kcat0: np.ndarray)

Biologically plausible bounds for proposed kcats, in 1/s.

Proposals outside these bounds are not hypotheses worth spending an FBA solve on. The window is 1e-2 to 1e4 for an ordinary kcat, and always widened far enough to contain the prior with a hundred-fold margin either side, so an unusually slow or fast enzyme keeps a window around itself rather than being clipped towards the generic range.

Parameters

kcat0: Prior kcat per tunable row, 1/s.

Returns

tuple of numpy.ndarray: Lower and upper bound per row.

load_evotune_data#

load_evotune_data(adapter: 'ModelAdapter')

Load the three evotune-kcat-tuning experimental data files.

Each file is optional and resolved as <adapter.params.path>/data/<filename>; a missing file yields None (flux_data/max_grate) or an empty list (zero_flux) rather than raising.

Parameters

adapter: A ModelAdapter whose params.path locates the data/ folder.

Returns

EvotuneData :

MaskSizePoint#

Tuned fit at one candidate target_impact_share.

Attributes

target_impact_share: The share this point was tuned at.
n_selected: Free parameters .tuning.select_tunable_mask chose at this share (tied isozyme copies count once).
rmse: The lowest value reached in each seed's EvotuneResult.rmse_trace -- not rmse_trace[-1], because rmse_trace only records plain RMSE at a step that improved the (possibly prior-penalised) objective, per cmaes_kcat_tuning's own docstring; with the default nonzero prior_penalty_weight a later generation can improve the objective while its plain RMSE is slightly worse than an earlier generation's, making the final entry not necessarily the best one recorded. min(rmse_trace) is the best that was recorded either way, though a candidate that improved plain RMSE without improving the objective is invisible to both and to rmse_trace itself -- set params.prior_penalty_weight=0 before a sweep meant to measure fit alone, if that residual gap matters. Kept per-seed rather than pre-averaged so a caller can look at the spread directly.
seeds: The seeds each entry of rmse came from.

movement_in_sigma#

movement_in_sigma(kcat: np.ndarray, kcat0: np.ndarray, sigma0_log: np.ndarray)

Per-kcat displacement from the prior, in prior standard deviations.

|log(kcat / kcat0)| / sigma0_log, matching MATLAB's devFromPrior = abs(muLog - log(kcat0)) ./ sigma0log.

n_changed#

n_changed(kcat: np.ndarray, kcat0: np.ndarray, rel_tol: float=0.02)

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

The 2% default matches the threshold MATLAB reports its per-source "unchanged" percentages at.

n_free#

n_free(tie_map: np.ndarray)

How many parameters remain free once ties are applied.

parsimony_frontier#

parsimony_frontier(score: Callable[[np.ndarray], float], kcat: np.ndarray, kcat0: np.ndarray, sigma0_log: np.ndarray, thresholds: Sequence[float]=(0.0, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0))

Score kcat with progressively more of it reverted to prior.

Parameters

score: Scores one kcat vector; one call per threshold.
kcat, kcat0, sigma0_log: Tuned vector, prior vector, and per-kcat prior sigma.
thresholds: Sigma cutoffs to sweep, ascending. 0.0 reverts nothing and so reproduces kcat's own score.

Returns

list of FrontierPoint, in the order given. :

revert_below#

revert_below(kcat: np.ndarray, kcat0: np.ndarray, sigma0_log: np.ndarray, threshold: float)

Put every kcat that moved less than threshold sigma back on

its prior, leaving the rest at their tuned values.

screen_kcat_leverage#

screen_kcat_leverage(model: 'EcModel', *, adapter: Optional['ModelAdapter']=None, params: Optional['EvotuneParams']=None, evotune_data: Optional[EvotuneData]=None, okp_method: Optional[str]=None, bio_rxn: Optional[str]=None, make_anaerobic: Optional[Callable[['EcModel'], None]]=None, change_protein_biomass: Optional[Callable[['EcModel', float], None]]=None, tunable_mask: Optional[np.ndarray]=None, fold: float=2.0, n_proc: Optional[int]=None)

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.tie_isozymes) is perturbed up and down by fold, and its leverage is the largest resulting change in plain RMSE -- always the plain fit, never the prior_penalty_weight-penalised objective, since leverage asks what the data supports, independent of how much a move is charged for. Groups are ranked by that leverage weighted by sigma0_log -- .parsimony.identifiability_mask's convention, where a trusted source needs proportionally more measured effect to rank alongside an untrusted one. This trust-weighting means the reported cum_leverage_share is a different quantity from .parsimony.impact_share, which is unweighted.

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. FVA-unreachable kcats are not filtered out beforehand: a kcat that cannot carry flux under any condition shows zero leverage here by construction and simply sorts to the bottom, which is a more direct test than a separately computed static reachability mask.

Returns

pandas.DataFrame: One row per tie-group, sorted by rank descending: rxn_id (the group's representative), n_isozymes, source_group, kcat0, leverage, sigma0_log, rank (leverage * sigma0_log), and cum_leverage_share (the running total of rank, as a fraction of its sum over every row). An internal _positions column (a tuple of indices into model.ec.rxns) is what select_tunable_mask needs to expand a row back into a mask; drop it before showing the table to a person.

select_tunable_mask#

select_tunable_mask(model: 'EcModel', screen: pd.DataFrame, *, target_impact_share: float=0.9)

Build a tunable_mask from a screen_kcat_leverage report.

Keeps the fewest highest-ranked groups whose combined leverage reaches target_impact_share of the total -- a relative cutoff, so the same target_impact_share 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. target_impact_share=0.9 is a starting point, not a validated universal constant; check the resulting mask's size against screen_kcat_leverage's own curve before trusting it on a new model.

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

sensitivity_tuning#

sensitivity_tuning(model: 'EcModel', *, desired_growth_rate: Optional[float]=None, fold_change: float=10.0, prot_to_ignore: Optional[list[str]]=None, max_iterations: int=1000, bio_rxn: Optional[str]=None, verbose: bool=True)

Iteratively bump the most-limiting kcat by fold_change until

the model can reach desired_growth_rate.

Each iteration:

Solves the LP at the current kcat values.
If growth didn't increase since the last iteration, warns and breaks (likely a substrate / pool constraint, not a kcat one).
If growth reached the target, returns.
Otherwise picks the usage_prot_* reaction with the highest flux (= most-produced enzyme), excluding any whose UniProt ID is in prot_to_ignore.
Finds the catalysed reaction consuming the most of that enzyme's pseudometabolite and multiplies that reaction's ec.kcat by fold_change.
Annotates ec.notes with the original kcat and source (unless the source is already "sensitivityTuning").
Re-applies kcat constraints for that single reaction.

Ported from GECKO MATLAB: src/geckomat/kcat_sensitivity_analysis/sensitivityTuning.m.

Parameters

model: EcModel with the protein pool / usage rxn machinery installed and a populated ec.kcat. Mutated in place.
desired_growth_rate: Target growth rate. Defaults to model.adapter.params.gr_exp.
fold_change: Per-iteration kcat multiplier.
prot_to_ignore: UniProt IDs whose usage rxns are excluded from being picked as the limiting enzyme.
max_iterations: Hard cap on tuning iterations. A backstop against a target growth rate that is unreachable by raising kcats alone.
bio_rxn: Biomass reaction id to set as the objective. Defaults to params.bio_rxn from the adapter when one is attached; if neither is available the model's current objective is used.
verbose: Whether per-iteration progress is logged at INFO.

Returns

TunedKcatsResult :

Raises

ValueError: If the initial LP is infeasible, or if desired_growth_rate is needed but no adapter / value is available.
NotImplementedError: If model.ec.gecko_light is True.

SigmaFitterResult#

Outcome of a sigma scan.

Attributes

sigma: The sigma value that minimised the absolute relative error between predicted and experimental growth.
sigma_grid: Sigma values that were tried.
growth_grid: LP-optimal growth at each sigma.
error_grid: Absolute relative error at each sigma, in percent.

simulate_evotune_dataset#

simulate_evotune_dataset(model: 'EcModel', flux_data: 'FluxData', *, constrain: bool, zero_flux_rxns: list[str], bio_rxn_id: str, make_anaerobic: Optional[Callable[['EcModel'], None]]=None, change_protein_biomass: Optional[Callable[['EcModel', float], None]]=None)

Simulate every condition in one tuning dataset.

Mirrors abc_max.m's rmsecal FBA half:

Every exchange reaction that appears as a condition name anywhere in flux_data.conds (matched by name against flux_data.exch_mets) is blocked (lower_bound = 0) up front, across the board -- not just the reaction for the row being solved -- so one row's carbon source can't leak flux into another row's simulation.
For row i, that row's own carbon-source reaction is unblocked: fixed at the measured uptake rate if constrain is True (evotune_data.flux_data), or fully opened to -1000 if constrain is False (evotune_data.max_grate -- "what's the best possible growth on this carbon source", not "growth at this measured rate").
Optional per-condition adjustments run if the corresponding callable is supplied: anaerobic switch when the oxygen exchange column reads exactly 0, protein-content rescaling when Ptot is given and nonzero.
Solve; record growth, plus (constrain datasets only, since max_grate has no exchange measurements to compare against) every non-NaN measured exchange's simulated flux and every zero_flux_rxns entry's simulated flux.

All mutation happens inside with model:, which cobra reverts on exit -- only bounds are touched here, so this is safe to call repeatedly against one persistent (per-worker) model instance; a per-particle EcModel.copy() is not used (prohibitively expensive on a real-scale model).

Parameters

model: EcModel with the current particle's kcats already applied (via apply_kcat_constraints).
flux_data: The dataset to simulate.
constrain: True for flux data (fix uptake at the measured rate), False for max-growth data (open uptake fully).
zero_flux_rxns: Reaction IDs assumed zero-flux in every condition (evotune_data.zero_flux); the row's own carbon source is skipped even if it's also listed here.
bio_rxn_id: The biomass reaction ID (adapter.params.bio_rxn).
make_anaerobic, change_protein_biomass: Optional adapter-specific hooks (MATLAB: modelAdapter.makeModelAnaerobic / .changeProteinBiomass). geckopy has no generic, organism-agnostic port of these yet -- in MATLAB they're written per model-adapter subclass (e.g. YeastGEMAdapter.m). Left as an injection point rather than a hard requirement: data with no anaerobic/Ptot rows never needs them, and passing None for either just skips that adjustment instead of raising.

Returns

list[ConditionSimResult]: One entry per row of flux_data.conds, in order.

Raises

ValueError: If any flux_data.conds entry can't be matched by name to an exchange metabolite in flux_data.exch_mets.

source_movement#

source_movement(kcat: np.ndarray, kcat0: np.ndarray, groups: np.ndarray, order: Sequence[str], near_tol: float=0.1)

Median fold change and near-prior fraction, per source group.

Two medians are reported per source, separating how many kcats a result moves from how far it moves them. median_fold is taken over every kcat of the source; where most of a source sits at its prior -- which a restrictive tunable_mask guarantees -- it saturates at 1.0 and carries no information about the tier. median_fold_moved conditions on the kcats that did move, and so stays informative at any sparsity.

Parameters

order: Source names from most to least trusted. Used only to report whether movement respects that ranking.
near_tol: A kcat within 1 + near_tol fold of its prior counts as untouched, and is excluded from median_fold_moved.

Returns

dict keyed by source name, plus ``"_ordered"``: whether :

```median_fold_movedis non-decreasing alongorder``, i.e.` :

whether the result moved trusted kcats less far than untrusted :

`ones. A source with nothing moved reports ``median_fold_moved``` :

1.0 and takes no part in that check. :

sweep_tunable_mask_size#

sweep_tunable_mask_size(model: 'EcModel', *, adapter: Optional['ModelAdapter']=None, params: Optional['EvotuneParams']=None, evotune_data: Optional['EvotuneData']=None, okp_method: Optional[str]=None, bio_rxn: Optional[str]=None, make_anaerobic: Optional[Callable[['EcModel'], None]]=None, change_protein_biomass: Optional[Callable[['EcModel', float], None]]=None, screen: Optional[pd.DataFrame]=None, target_impact_shares: Sequence[float]=(0.3, 0.5, 0.7, 0.8, 0.9, 0.95), seeds: Sequence[int]=(0, 1), n_proc: Optional[int]=None, verbose: bool=True)

Tune at each of target_impact_shares and report the fit.

The screen is computed once (or reused, via screen) and shared across every share and seed; only .tuning.cmaes_kcat_tuning reruns per point. Each run starts from model's own current ec.kcat -- model is not mutated (a fresh copy is tuned each time), unlike .tuning.cmaes_kcat_tuning itself.

Cost is one screen (comparable to a full tuning run) plus len(target_impact_shares) * len(seeds) tuning runs, each sized by however many parameters that share selects -- the later points in an increasing sequence of shares cost more, both because CMA-ES scales its population with dimension and because a wider search space usually takes longer to explore.

The defaults are a coarse scout, not a guarantee of finding the true optimum: five or six points cannot resolve a curve whose real knee sits inside one of the gaps between them, and how wide that gap needs to be to matter is exactly what a first sweep is for finding out. Look at n_selected after a first sweep -- if consecutive shares jump by a large multiple in parameter count (target_impact_share is closer to the leverage tail than to its bulk there, e.g. the jump from 0.9 to 0.95 was 5x on one real model), or if the recommendation lands next to the widest such gap, that gap is where the true knee could be hiding; rerun with a few shares added inside it. The same applies at the low end: if the smallest share measured isn't clearly worse than the best one, it has not necessarily bracketed the true minimum, and a smaller share is worth adding to confirm it.

seeds defaults to two, not one, doubling the cost of the default call: recommend_target_impact_share's noise-floor check needs at least two seeds at a point to have anything to check against, and with only one it silently falls back to a cruder rule everywhere. Pass a single seed for the very first, cheapest possible scout of the shape; widen to more once a promising region is known, since that is where seed noise is most likely to matter (fit values close enough together for it to decide the comparison).

Parameters

model: EcModel with a populated ec.kcat, used as the tuning starting point at every share/seed. Not mutated.
target_impact_shares: Candidate shares to measure, any order (sorted on return).
seeds: Seeds tuned at every share.
screen: A precomputed .tuning.screen_kcat_leverage report. Computed once, fresh, if not given.

adapter, params, evotune_data, okp_method, bio_rxn, make_anaerobic, :

change_protein_biomass, n_proc, verbose: Forwarded to .tuning.screen_kcat_leverage and .tuning.cmaes_kcat_tuning; see their docstrings.

Returns

list of :class:MaskSizePoint, one per share, sorted by :

```target_impact_share`` ascending.` :

truncate_values#

truncate_values(arr: np.ndarray, columns: Optional[list[int]]=None)

Round values in arr to ~6 significant figures.

For each value v:

order_magn = max(ceil(log10(|v|)), 0)
result     = round(v, 6 - order_magn)

Zero, NaN, and infinite values are returned unchanged.

For 2D input, only the columns selected by columns are processed (default: all). Other columns are copied through unchanged. The input is never mutated; a new array is returned.

Ported from GECKO MATLAB: src/geckomat/kcat_sensitivity_analysis/truncateValues.m.

Parameters

arr: 1-D or 2-D numeric array.
columns: Column indices to process for 2-D input. Default: every column. Ignored for 1-D input.

Returns

numpy.ndarray: New array with the truncated values.

tune_prior_penalty_weight#

tune_prior_penalty_weight(model: 'EcModel', *, adapter: Optional['ModelAdapter']=None, params: Optional['EvotuneParams']=None, evotune_data: Optional[EvotuneData]=None, okp_method: Optional[str]=None, bio_rxn: Optional[str]=None, make_anaerobic: Optional[Callable[['EcModel'], None]]=None, change_protein_biomass: Optional[Callable[['EcModel', float], None]]=None, tunable_mask: Optional[np.ndarray]=None, screen: Optional[pd.DataFrame]=None, target_impact_share: float=0.9, candidates: Sequence[float]=(0.0, 0.01, 0.03, 0.1), seeds: Sequence[int]=(0, 1), fold_threshold: float=2.0, popsize: Optional[int]=None, n_proc: Optional[int]=None, verbose: bool=True)

Sweep prior_penalty_weight and report the fit/reproducibility

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

prior_penalty_weight 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, giving up fit for nothing. This runs cmaes_kcat_tuning at every value in candidates, each at every seed in seeds, all against the same tunable set (built once, the way cmaes_kcat_tuning would on its own if tunable_mask/screen aren't given), and reports both fit and cross-seed reproducibility per value.

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

model is restored to its original kcats before every run and left that way when this returns -- unlike cmaes_kcat_tuning, which mutates model with its result, this function only reports, since there is no single "best" candidate for it to leave the model in. Fewer than two seeds makes every reproducibility column NaN, since there is nothing to compare.

Parameters

candidates: prior_penalty_weight values to try.
seeds: Seeds to run each candidate at. Reproducibility columns pool every pair of seeds, so more than two sharpens them but costs proportionally more.
fold_threshold: A kcat counts as "moved" past this fold change from its prior. Matches .parsimony.fold_change's convention.

tunable_mask, screen, target_impact_share, popsize, n_proc, :

make_anaerobic, change_protein_biomass, okp_method, bio_rxn: As in cmaes_kcat_tuning.

Returns

pandas.DataFrame: One row per candidate: prior_penalty_weight, n_seeds, distance_mean and distance_spread (max minus min across seeds), fit_cost (distance_mean relative to the best candidate's, as a fraction), n_changed_mean, n_movers (kcats past fold_threshold in either seed of a pair, summed over every pair), n_both_moved (past it in both -- the denominator for direction agreement, since a kcat only one seed moved has no direction to agree about), pct_direction_agree, median_fold_spread and max_fold_spread (among movers, the largest fold-change between what the two seeds landed on).

TunedKcatsResult#

Outcome of a sensitivity-tuning run.

Attributes

rxns: IDs of the reactions whose kcat was bumped, in tuning order.
rxn_names: Reaction names parallel to rxns (falls back to the reaction ID when it has no name).
enzymes: For each tuned reaction, the catalysing enzyme ID(s), ;-joined when more than one enzyme is involved.
old_kcat: kcat value before tuning, parallel to rxns.
new_kcat: kcat value after tuning, parallel to rxns.
source: Original kcat provenance (model.ec.source) before tuning, parallel to rxns.

kcat sensitivity tuning (Python)#

Functions#

Reference#

annotate_from_model#

apply_ties#

base_reaction#

best_parsimonious#

build_sigma0_log#

classify_kcat_source#

classify_kcat_sources#

cmaes_kcat_tuning#

compute_excarbon#

ConditionSimResult#

Correction#

corrections#

corrections_tsv#

dataset_rmse#

evotune_distance#

EvotuneData#

EvotuneResult#

find_max_value#

fit_sigma#

fold_change#

FrontierPoint#

identifiability_mask#

impact_share#

isozyme_tie_map#

kcat_bounds#

load_evotune_data#

MaskSizePoint#

movement_in_sigma#

n_changed#

n_free#

parsimony_frontier#

recommend_target_impact_share#

revert_below#

screen_kcat_leverage#

select_tunable_mask#

sensitivity_tuning#

SigmaFitterResult#

simulate_evotune_dataset#

source_movement#

sweep_tunable_mask_size#

truncate_values#

tune_prior_penalty_weight#

TunedKcatsResult#