API Reference

fit_rate_equation(fp::FittingProblem, optimizer;
    n_restarts=20, maxtime=60.0, maxiters=10_000_000,
    abstol=nothing, reltol=nothing, callback=nothing,
    lb=fill(-15.0, length(fitted_params(fp.mechanism))),
    ub=fill(15.0, length(fitted_params(fp.mechanism))),
    solver_kwargs=(;))

Fit rate constants by minimizing loss! using Optimization.jl.

Runs n_restarts independent optimizations from random initial points and returns the best result.

maxtime and maxiters are Optimization.jl common solver options forwarded to every Optimization.solve call. abstol, reltol, and callback are common options forwarded only when set (otherwise each solver uses its own default). solver_kwargs is a NamedTuple of solver-specific options forwarded verbatim to solve (e.g. (; popsize=200) for a CMA-ES solver that supports it). A key present in both a named common option and solver_kwargs takes the solver_kwargs value. lb/ub are the log-space bounds passed to the OptimizationProblem.

Rescaling is driven by fp.scale_k_to_kcat: when it is a Real, the returned parameters are rescaled so that _kcat_forward(mechanism, params) ≈ fp.scale_k_to_kcat; when it is nothing, the raw (unrescaled) parameters are returned (the data fixes the absolute scale).

Returns a NamedTuple (params, loss, retcode) where:

• params: fitted rate constants as a NamedTuple
• loss: the best loss value achieved
• retcode: the Symbol form of the best restart's sol.retcode. Only :Success indicates the optimizer converged on its own criteria; any other value (e.g. :MaxTime — hit the time budget; :Failure; :Default — ran but did not flag success) means the fit should be treated as un-converged (check retcode !== :Success).

identify_rate_equation(prob; optimizer,
    min_beam_width=50, loss_rel_threshold=1.3, loss_abs_threshold=0.001,
    loss_parsimony_threshold=0.99,
    max_param_count=20, n_restarts=20, maxtime=60.0, maxiters=10_000_000,
    abstol=nothing, reltol=nothing, callback=nothing, solver_kwargs=(;),
    n_cv_candidates=5, se_threshold=1.0, perm_p_threshold=0.16,
    save_dir=_default_save_dir(), show_progress=true)

Find the best rate equation for the given reaction and data using beam search.

Keyword Arguments

• min_beam_width::Int = 50: minimum mechanisms to keep per param-count tier
• loss_rel_threshold::Float64 = 1.3: relative tolerance for beam selection (see "Beam selection" below)
• loss_abs_threshold::Float64 = 0.001: absolute tolerance for beam selection
• loss_parsimony_threshold::Float64 = 0.99: a mechanism keeps expanding only if its loss is within this factor of the best model of any smaller parameter count — an added parameter must earn its keep. Combined with the other loss thresholds via min; min_beam_width is a cumulative per-count budget (see "Beam selection" below). Inf disables it.
• max_param_count::Int = 20: stop expanding beyond
• eq_complexity_filter::Int = 337: skip (before fitting, before any derivation) any mechanism whose rate equation is more complex than this — roughly the number of terms in the rate equation's denominator: V×τ, the number of RE-segments times the spanning-tree count of the catalytic segment graph, i.e. the products the compiled equation evaluates on every call. The default 337 is one above a fully steady-state random-order bi-bi (V×τ = 336), so such a mechanism passes and anything more complex is skipped — equations denser than that are impractical to fit and can blow up derivation/codegen. Computed from the mechanism graph alone.
• optional_allosteric_regulators::Vector{Symbol} = Symbol[]: declared allosteric regulators the seed may omit. A reaction that declares regulators seeds the beam fully regulated by default; naming a regulator here lets the beam also start from seeds that do not bind it (added back later by expansion). Listing every declared allosteric regulator restores the unregulated init_mechanisms seed.
• optional_competitive_inhibitors::Vector{Symbol} = Symbol[]: declared competitive inhibitors the seed may omit, as above.
• optimizer: Optimization.jl optimizer (required). Recommended: CMAEvolutionStrategyOpt() from OptimizationCMAEvolutionStrategy.
• n_restarts::Int = 20: multi-start restarts per fit
• maxtime::Real = 600.0: max time per fit (seconds; common solver option, forwarded to Optimization.solve)
• maxiters::Integer = 10_000_000: max iterations per optimizer run (common solver option, forwarded to Optimization.solve)
• abstol/reltol/callback = nothing: Optimization.jl common solver options, forwarded to Optimization.solve only when set
• solver_kwargs::NamedTuple = (;): solver-specific options forwarded verbatim to Optimization.solve (e.g. (; popsize=200) for a CMA-ES solver that supports it); the caller matches its contents to optimizer
• n_cv_candidates::Int = 5: LOOCV top N unique-rate-equation candidates per param count
• se_threshold::Float64 = 1.0: paired 1-SE multiplier for model selection. Simpler-model bucket accepted iff its mean paired loss difference vs the best bucket is ≤ se_threshold * std(diffs)/sqrt(n_folds). Default 1.0 is the textbook "1-SE rule".
• perm_p_threshold::Float64 = 0.16: minimum one-sided permutation p-value for model selection. Simpler-model bucket accepted iff p > perm_p_threshold under the sign-flip null. Default 0.16 matches paired 1-SE empirically.
• save_dir::String = _default_save_dir(): output directory for the search CSVs (initial_mechanisms.csv + equation_search_iteration_N.csv), plus loocv_results.csv (the LOOCV table for every cross-validated candidate) and best_equation.csv (the selected equation and its fitted parameters)

Beam selection

A mechanism at parameter count n qualifies for the next-level beam if either:

• its loss ≤ min(loss_rel_threshold * best(n) + loss_abs_threshold, loss_parsimony_threshold * best(<n)), where best(n) is the lowest loss at parameter count n and best(<n) is the lowest loss over all smaller counts; the second term is dropped at the base count (no smaller level exists yet),
• OR it falls under the min_beam_width budget: a cumulative per-count allowance that expands at least min_beam_width mechanisms at each count over the whole search, then stops — spent once, not re-granted each time the count is revisited.

The additive term protects against best_loss approaching zero (simulated / very-low-loss data) where a purely multiplicative threshold would collapse the beam to the single best mechanism.

Model selection (LOOCV)

For each n_params bucket below n_min (lowest mean fold-loss) the rule computes paired fold-loss differences between the bucket's representative and n_min's, then accepts the simpler bucket iff BOTH:

1. Paired 1-SE rule: mean(diffs) ≤ se_threshold * std(diffs)/sqrt(n_folds).
2. One-sided permutation test: permutation_p > perm_p_threshold, where p is computed by exact enumeration when n_folds ≤ 20 and Monte Carlo (10⁶ samples) otherwise.

Returns the smallest passing n_params; falls through to n_min if none pass. Within the chosen bucket the mechanism with lowest training loss wins. Per-bucket representative = the row with the lowest cv_score in that bucket. Diagnostic columns mean_loss_diff, se_paired, permutation_p are surfaced in cv_results; the n_min bucket has 0.0 in all three.

metabolites(m::EnzymeMechanism) → Tuple{Symbol,...}

Return distinct metabolite names (substrates ∪ products ∪ regulators) as a tuple of Symbols in declaration order, deduplicated. The fitter uses this as the key set for the per-datapoint concentration NamedTuple passed to rate_equation.

julia> using EnzymeRates

julia> m = @enzyme_mechanism begin
           substrates: S
           products: P
           steps: begin
               E + S <--> E(S)
               E(S) <--> E + P
           end
       end;

julia> metabolites(m)
(:S, :P)

parameters(m::EnzymeMechanism, [mode])
parameters(m::AllostericEnzymeMechanism, [mode])

Return the parameter names required for the given mode as a tuple of Symbols.

Modes

• Reduced (default): independent k's + Keq + E_total. The set of symbols the user supplies to evaluate the Haldane-reduced rate equation. Returned for both EnzymeMechanism and AllostericEnzymeMechanism.
• Full: all raw rate-constant symbols + E_total. For EnzymeMechanism this is "all 2N k's + E_total." For AllostericEnzymeMechanism it composes the catalytic raw A-state symbols + every I-state mirror (catalytic + regulatory + synthesized dep) + reg-site A-state K's (skipping :OnlyI ligands) + :L + :E_total. The allosteric Full mode enumerates the complete raw-symbol set; no rate_equation method is defined for (::AllostericEnzymeMechanism, ::FullMode), so this mode is for symbol enumeration, not runtime evaluation.

rate_equation(m, concs, params, [mode])

Return the net reaction rate of mechanism m at the metabolite concentrations concs and parameters params (both NamedTuples). The rate equation itself is derived from the mechanism's elementary steps by the King–Altman/Cha method: rapid-equilibrium steps collapse to binding constants, steady-state steps are assembled into the King–Altman determinant, and the two combine into the quasi-steady-state flux through the whole mechanism — so one call returns the overall turnover, not the rate of a single step.

mode selects how parameters enter the equation. The default Reduced applies the Haldane/Wegscheider reduction, deriving the dependent rate constants from Keq and the independent parameters; Full instead takes every rate constant as independent.

The derivation runs once, at compile time: the body is generated as a single arithmetic expression with no allocations, loops, or matrix operations. Use rate_equation_string to inspect the derived equation symbolically.

rate_equation_string(m, [mode]) -> String

Return the symbolic rate equation for mechanism m as a multi-line String (it returns the text — it does not print). mode is Reduced (default) or Full; pass a concrete Mechanism / AllostericMechanism or its compiled EnzymeMechanism singleton.

The string is a runnable transcript of how rate_equation evaluates: a (; …) = params destructure line, a (; …) = concs destructure line, then the v = E_total * (num) / (den) line. In Reduced mode, dependent rate constants are listed first under # Wegscheider constraints: and # Haldane constraints: headers — the thermodynamic identities that eliminate parameters — and only the independent set appears in the params destructure. In Full mode every rate constant is independent, so there is no constraint section. Full mode is defined for EnzymeMechanism only; an AllostericEnzymeMechanism supports Reduced mode only.

Use print on the result to see the multi-line layout without escaped newlines.

julia> using EnzymeRates

julia> m = @enzyme_mechanism begin
           substrates: S
           products: P
           steps: begin
               E + S ⇌ E(S)
               E(S) <--> E(P)
               E(P) ⇌ E + P
           end
       end;

julia> print(rate_equation_string(m))
(; K_P_E, K_S_E, k_ES_to_EP, Keq, E_total) = params
(; S, P) = concs
# Haldane constraints:
k_EP_to_ES = (1 / Keq) * K_P_E * (1 / K_S_E) * k_ES_to_EP
v = E_total * (k_ES_to_EP * S / K_S_E - k_EP_to_ES * P / K_P_E) / (1 + P / K_P_E + S / K_S_E)

rescale_parameter_values(m, params::NamedTuple; scale_k_to_kcat=1.0)

Rescale SS rate constants so that _kcat_forward(m, result) ≈ scale_k_to_kcat. Non-SS parameters (K's, Keq, E_total, L, regulatory K's) are unchanged.

@allosteric_mechanism begin
    substrates: F6P
    products:   F16BP
    catalytic_multiplicity: 2
    allosteric_regulators: A::OnlyA, I::OnlyI

    catalytic_steps: begin
        E + F6P ⇌ E(F6P)        :: EqualAI
        E(F6P) <--> E(F16BP)    :: EqualAI
        (E(F16BP) ⇌ E + F16BP)  :: EqualAI
    end

    regulatory_site(multiplicity = 4): begin
        ligands: A
    end
    regulatory_site(multiplicity = 4): begin
        ligands: I
    end
end

Build an AllostericEnzymeMechanism (MWC, two conformations).

• substrates:, products:, catalytic_inhibitors: accept comma-separated bare symbols.
• allosteric_regulators: requires name::Tag per entry, where Tag is one of OnlyA, OnlyI, EqualAI, NonequalAI.
• catalytic_multiplicity: N is the subunit count for the catalytic site (default 1).
• catalytic_steps: begin ... end is required (exactly once); each step or parenthesized step-group must carry a ::Tag from the same set. Function- call species notation (E(F6P), Estar(B; residual = A - P)) is supported.
• regulatory_site(multiplicity = N): begin ligands: L1, L2 end declares one regulatory site per block with multiplicity N and the ligands listed inside. Each ligand must appear in allosteric_regulators:; a ligand may not appear in two sites. Ligands declared in allosteric_regulators: but not assigned to any regulatory_site(...): block default to a single-ligand site at the catalytic multiplicity.

@enzyme_mechanism begin
    substrates: S
    products:   P
    regulators: I

    steps: begin
        E + S ⇌ E(S)                   # function-call species notation
        (E(S) ⇌ E(P), E_alt(S) ⇌ E_alt(P))   # parenthesized → shared kinetics
        E(S) + I ⇌ E(S, I)             # dead-end
        E(P) ⇌ E + P
    end
end

Build a plain (non-allosteric) EnzymeMechanism.

• substrates:, products:, regulators: accept comma-separated bare symbols. Atom brackets (S[C]) are rejected at the mechanism level.
• regulators: entries are treated as CompetitiveInhibitors when later passed to EnzymeReaction.
• Same-kinetics groups are expressed via parenthesized step-groups; no constraints: block needed.
• Allosteric-only constructs (site(...) / ::Tag / allosteric_regulators: / catalytic_inhibitors:) are rejected.

Species notation on step sides:

• Bare Symbol that matches a declared metabolite → that metabolite.
• Bare Symbol otherwise (e.g. E, Estar, ES) → conformation-only species named after the Symbol.
• E(S) / E(S, P) → species with conformation :E and bound metabolites; synthesized name is :E_<bound...> with bound names sorted alphabetically (matching name(::Species)).
• Estar(; residual = A - P) → species with empty bound and a residual recording +A / −P; synthesized name is :Estar_res_+A_-P.
• Estar(B; residual = A - P) → bound + residual; name :Estar_B_res_+A_-P.

Conformation labels cannot shadow declared metabolite names.

@enzyme_reaction begin
    substrates: S[C6H12O6], ATP[C10H16N5O13P3]
    products:   P[C6H13O9P]
    competitive_inhibitors: I            # bare OK; mults default to catalytic
    allosteric_regulators: A(1, 2, 4)    # per-reg mults required
    allowed_catalytic_multiplicities: (1, 2, 4)
    # or shorthand:
    # oligomeric_state: 2
end

Emit an EnzymeReaction.

• substrates: / products: — comma-separated entries with required atom brackets (S[C6H12O6]). Multi-atom forms like [C2,N] are allowed.
• competitive_inhibitors: / dead_end_inhibitors: — CompetitiveInhibitor entries (catalytic-site binding). May be bare I (multiplicities default to allowed_catalytic_multiplicities) or I(m1, m2, ...) to override.
• allosteric_regulators: — AllostericRegulator entries. Each entry must declare per-regulator multiplicities as A(m1, m2, ...) or a single value A(m). An entry may carry a type tag, A::Activator or A::Inhibitor; untagged entries default to :unspecified. Type tags are only valid on allosteric_regulators: entries.
• allowed_catalytic_multiplicities: — tuple of positive Ints (default (1,)).
• oligomeric_state: N — shorthand for allowed_catalytic_multiplicities: (N,).
• shared_catalytic_site: (sub, prod), ... — pairs of substrate/product names (in either order) that bind at the same catalytic site.

AllostericEnzymeMechanism(am::AllostericMechanism)

Lift an AllostericMechanism to its singleton AllostericEnzymeMechanism type. The catalytic side becomes an EnzymeMechanism lifting through Mechanism(am.reaction, am.cat_steps). Catalytic and regulatory allosteric data are encoded directly into the type parameters.

EnzymeMechanism{Sig}

Singleton type encoding an enzyme mechanism. Sig is the tuple (reaction_sig, steps_sig) produced by _sig_of(::Mechanism):

• reaction_sig encodes the EnzymeReaction (reactants, regulators, catalytic multiplicities) as a tuple of tuples of Symbols/Ints.
• steps_sig encodes Vector{Vector{Step}} as a nested tuple where the outer level is kinetic groups and the inner level is Step data (from-species, to-species, bound metabolite, is-equilibrium).

The conversion functions are _sig_of and _mechanism_from_sig; the exact layout is internal — users construct via EnzymeMechanism(::Mechanism).

EnzymeReaction

The public concrete reaction descriptor: substrates, products, optional regulators, and the set of catalytic multiplicities the mechanism enumerator is allowed to consider. It is the entry point to the package — pair it with rate data in an IdentifyRateEquationProblem, or pass it to EnzymeRates.init_mechanisms to enumerate candidate mechanisms.

Fields

• reactants::Vector{ReactantAtoms} — substrates and products, each carrying its per-atom inventory (used for ping-pong residual bookkeeping and atom conservation across steps).
• regulators::Vector{RegulatorMults} — competitive inhibitors and allosteric regulators, each with its allowed oligomeric multiplicities.
• allowed_catalytic_multiplicities::Vector{Int} — the oligomeric states the allosteric enumerator may assign to the catalytic core.
• shared_catalytic_site::Vector{Tuple{Symbol,Symbol}} — (substrate, product) pairs that bind the same catalytic site, and so are never both bound to it at once.

Construct one with the @enzyme_reaction DSL. Each reactant takes an atom bracket (S[C], A[C1H1]); the brackets are load-bearing for ping-pong and multi-substrate reactions. Reactants and regulators are sorted by name in the constructor, so two equivalent declarations compare equal under ==/hash.

julia> using EnzymeRates

julia> rxn = @enzyme_reaction begin
           substrates: S[C]
           products:   P[C]
       end;

julia> rxn isa EnzymeReaction
true

julia> rxn
EnzymeReaction: S ⇌ P

FittingProblem{M, D}

Holds pre-processed experimental data and mechanism info for fitting rate constants.

Fields

• mechanism: an AbstractEnzymeMechanism instance
• data: NamedTuple of column vectors (via Tables.columntable)
• group_point_indexes: row indices grouped by unique group values
• Keq: fixed equilibrium constant
• scale_k_to_kcat: a positive Real selects relative mode (per-group-centered loss); nothing selects absolute per-enzyme-turnover mode (uncentered loss)
• log_abs_rates: pre-computed log.(abs.(Rate))
• log_ratios_buffer: pre-allocated working buffer for loss computation

FittingProblem(mechanism::AbstractEnzymeMechanism, table;
    Keq::Real, scale_k_to_kcat::Union{Real,Nothing}=1.0)

Construct a FittingProblem from an enzyme mechanism and tabular data.

The table must have columns: group, Rate, and one column per metabolite matching metabolites(mechanism). Uses Tables.columntable for conversion.

scale_k_to_kcat selects the loss mode: a positive Real (default 1.0) treats the data as relative (per-group-centered loss); nothing treats it as absolute per-enzyme turnover (uncentered loss).

Rate values must be nonzero (zero rates produce -Inf in log space).

IdentifyRateEquationProblem{R, D}

Holds the reaction, experimental data, and equilibrium constant for rate equation identification.

Fields

• reaction: the EnzymeReaction instance
• data: NamedTuple of column vectors with :group, :Rate, and metabolite columns
• Keq: fixed equilibrium constant
• scale_k_to_kcat: target kcat for SS-rate rescaling before fitting (nothing = no rescaling, positive Float64 = target)

IdentifyRateEquationResults

Results from identify_rate_equation.

Fields

• best: the best AbstractEnzymeMechanism (lowest loss at optimal param count)
• cv_results: DataFrame with LOOCV results for top candidates per param count

Singleton instance for full (raw) mode.

Singleton instance for reduced (Haldane-Wegscheider) mode.