Modeling 3: estimation and identifiability

A model is only as good as its parameters, and data rarely determine all of them. In this tutorial you ask, before measuring anything, which parameters of the batch culture from Your first model a biomass measurement can determine; draw synthetic data; estimate the parameters with their uncertainty; and keep the result as a parameter set the store pins to the model and the data. Everything runs on the host, in Rust.

Do Your first model first: this page uses models/monod_batch and its scenario batch.

1. Ask the design first

Suppose you can measure biomass once an hour with a standard deviation of 0.05 g/L. Which of mu_max, k_s and y_xs can that determine?

labinetix analyze identifiability models/monod_batch/scenarios/batch.toml \
  --params mu_max,k_s,y_xs --sigma biomass=0.05

parameter  delta_msqr
mu_max     9.4561e1
k_s        1.5326e0
y_xs       5.6697e1
collinearity 5.557, determinant 4.5786e1
identifiable subset (collinearity <= 10): mu_max, y_xs, k_s

The measures are those of Brun, Reichert & Künsch (2001), from the sensitivities of the observed outputs to each parameter, scaled by the parameter and by the measurement’s standard deviation:

• delta_msqr, the importance: how much the observations move when the parameter moves by its own size, in standard deviations. mu_max and y_xs dominate; k_s matters about sixty times less, because it only shapes the few hours before the substrate runs out.
• collinearity γ: whether some combination of the parameters leaves the observations unchanged. Above 10 to 20, the data cannot tell them apart, whatever the noise.
• determinant ρ: the joint information of the set, a single number to compare designs.
• identifiable subset: the parameters, taken in order of importance, whose collinearity stays below the limit (--gamma-max, 10 by default).

These values hold at the scenario’s parameters. Rerun with other values in the scenario’s [parameters], other outputs in --sigma, or a scenario that samples more densely around the end of the batch, and compare: that is designing the experiment.

2. Draw synthetic data

To test the estimation before real data exist, draw data from the model itself, with seeded noise, as a dataset entity of the store:

labinetix analyze synthesize models/monod_batch/scenarios/batch.toml \
  --sigma biomass=0.05 --seed 1 --dataset batch-synthetic
labinetix define sync dataset/batch-synthetic

This writes models/monod_batch/datasets/batch-synthetic/: the CSV under data/ (t,output,value,sigma, one observation per line) and an authoring that records the provenance and pins the CSV by its sha256. define sync locks it; define check refuses a CSV that no longer matches. Real data take the same shape: put the CSV under data/ and write the authoring with its hash. A dataset too large for git can leave the CSV out; the pin stays.

# Synthetic: the prediction of models/monod_batch/scenarios/batch.toml with Gaussian noise (biomass=0.05), seed 1.
t,output,value,sigma
0.0,biomass,0.040829808447552234,0.05
1.0,biomass,0.11834098489341496,0.05
2.0,biomass,0.22948691093706208,0.05
…

3. A scenario to estimate on

Estimation integrates the model and its sensitivities to every parameter many times, at parameter values far from the defaults. An explicit method like rk4 has a stability limit that depends on the parameters: near the end of the batch the substrate decays at a rate of mu_max · x / (k_s · y_xs), and at a small k_s that is fast enough to make the scenario’s RK4 substeps unstable. The state survives, held at zero by its bound, but the sensitivities do not, and the estimate refuses to report statistics computed from them:

labinetix: the sensitivities at the fit are not usable (standard errors [3.3997717618147165e-20, 3.0050020840821523e-21, 0.0]): an explicit solver step may be unstable at the fitted parameters; estimate on a scenario with method = "ros2", or with more substeps

ros2 is a stiff method (Verwer et al. 1999), stable at any rate. Copy the scenario to models/monod_batch/scenarios/batch-ros2.toml and change its solver:

[solver]
method = "ros2"
step = 0.1
substeps = 4
steps = 120
record_every = 10

rk4 and euler are the methods a target runs; ros2 is the host’s, for analysis.

4. The experiment

An experiment says which scenario, which dataset, and which parameters to estimate within which bounds. Write models/monod_batch/experiments/batch.toml:

[experiment]
scenario = "batch-ros2"
dataset = "batch-synthetic"
estimate = ["mu_max", "k_s", "y_xs"]
doc = "Growth rate, Monod constant and yield from hourly biomass samples."

[bounds]
mu_max = [0.05, 2.0]
k_s = [0.001, 5.0]
y_xs = [0.1, 1.0]

[fit]
starts = 20
seed = 1

[mcmc]
samples = 4000
burn_in = 1000
seed = 2

Parameters not listed keep the scenario’s values. An observation may be at any time within the scenario’s horizon: one between two steps adds an event step to it. The chain must keep at least ten samples per estimated parameter after the burn-in.

5. Estimate

labinetix analyze estimate models/monod_batch/experiments/batch.toml --write
labinetix define sync parameter_set/batch

Without --write the estimate is a dry run: it prints the parameter set and says whether it differs from the committed authoring, which it leaves alone. Run again later, with nothing changed, it reproduces the committed set:

dry run: parameter_set/batch is identical to the committed models/monod_batch/parameter_sets/batch/authoring/batch.toml; --write replaces it

The estimate takes about a minute and a half on a Raspberry Pi 5, most of it the MCMC.

What happens:

1. Weighted least squares in the logarithm of the parameters, so they stay positive and within their bounds: minimize χ² = Σ ((y − value) / σ)².
2. Levenberg–Marquardt from the scenario’s values and from 19 random points within the bounds, with the exact sensitivities of the outputs to the parameters.
3. The Fisher information at the best fit gives the standard errors and correlations (Cramér–Rao, σ known), and the identifiability of section 1, now at the fitted values.
4. MCMC, because [mcmc] is stated: adaptive Metropolis started at the fit with the Fisher covariance, a uniform prior within the bounds.

The result is the parameter set entity models/monod_batch/parameter_sets/batch/. Its authoring holds everything, and pins what it came from:

[model]
entity = "model_scaffold/monod_batch"
version = "0.1.0"
content_hash = "sha256:8cd6…"

[dataset]
entity = "dataset/batch-synthetic"
version = "1.0.0"
content_hash = "sha256:9642…"

[experiment]
name = "batch"
content_hash = "sha256:d59a…"

[scenario]
name = "batch-ros2"
content_hash = "sha256:029a…"

[values]
mu_max = 0.4950322083020979
k_s = 0.02283696800910036
y_xs = 0.4995393236464814

[fit]
estimated = ["mu_max", "k_s", "y_xs"]
standard_errors = [0.022215099784811166, 0.3688345673900873, 0.0024160035599216454]
correlation = [[1.0, 0.996176238984639, 0.3735182442468857], [0.996176238984639, 1.0, 0.3756343622501496], [0.3735182442468857, 0.3756343622501496, 0.9999999999999999]]
chi2 = 9.076072903934577
n = 13
aic = 15.076072903934577
bic = 16.77092097631919
starts = 20
hits = 14

[identifiability]
parameters = ["mu_max", "k_s", "y_xs"]
delta_msqr = [70.74308698579979, 0.1967460725013617, 61.879206719127204]
collinearity = 16.184768558146317
determinant = 15.216725479842935
gamma_max = 10.0
subset = ["mu_max", "y_xs"]

[posterior]
samples = 4000
burn_in = 1000
acceptance = 0.1445
mean = [0.4955970456578878, 0.021451926187563686, 0.49794983324305864]
sd = [0.0027466628359791912, 0.029944138000744684, 0.002837483239080114]
q025 = [0.4917806897230583, 0.0011377876526486104, 0.49192451250360236]
q975 = [0.5045987365129262, 0.12445792927292575, 0.5031650284714936]

6. Read the result

The data were drawn at mu_max = 0.5, k_s = 0.1, y_xs = 0.5.

• mu_max and y_xs are determined: 0.495 ± 0.022 and 0.4995 ± 0.0024, within a fraction of a standard error of the truth.
• k_s is not: 0.023, with a standard error of 0.37, sixteen times the value. Its correlation with mu_max is 0.996: a smaller k_s with a slightly smaller mu_max fits the hourly samples as well. At the fit the collinearity is 16, and k_s drops out of the identifiable subset. Section 1 found all three identifiable at k_s = 0.1: identifiability depends on where in parameter space you are.
• The posterior says the same without a linear approximation: the 95 % interval of k_s runs from the lower bound 0.001 to 0.124 and contains the truth; those of mu_max and y_xs are narrow. Where the Fisher information and the posterior disagree, as for k_s, the posterior is the one to believe: the problem is not linear there.
• hits = 14: 14 of 20 starts ended at the same optimum (χ² within 0.1). Few hits mean a rugged cost; raise starts.
• chi2 = 9.1 for n − k = 10 degrees of freedom is what a correct model with the right σ gives. A χ² far above n − k means the model or the noise level is wrong.
• AIC and BIC compare models fitted to the same data: lower is better, and the penalty for parameters keeps the comparison fair.
• A composed model’s set adds [derived]: every coefficient the elements determine (nu_growth_oxygen = { value = …, standard_error = … }) at the fit, its standard error propagated from the estimated yields. Only the given coefficients are estimated.

The fix for k_s is a design question, back to section 1: sample more often near the end of the batch, or measure the substrate.

7. The parameter set stays honest

The parameter set pins the model, the dataset, the experiment and the scenario. Change any of them and define check refuses the parameter set until you re-estimate. Change the experiment’s [fit] seed to 3 and check:

labinetix define check

labinetix: define check: invalid entity state: parameter_set/batch: …/models/monod_batch/experiments/batch.toml is sha256:aa05136844992e3a548a9f046809ac353fefb4c061626607e2154a102c27724b, and it pins sha256:d59aa9ab21903ac7a01441b8500e7548562535b48a04c1ee3bac4d2ff0fd306b: re-run `labinetix analyze estimate models/monod_batch/experiments/batch.toml --write` and sync it

Set the seed back to 1 and define check passes again.

Re-estimating with --write rewrites the authoring, keeping its version; raise contract_version by hand when you publish a new set, and labinetix define tag parameter_set/batch tags it like any entity. Committed parameter sets reproduce byte for byte on the same platform: the test suite re-runs every experiment.

8. Several runs

Real data come in runs: cultivations with their own feed, their own start and their own samples. A dataset of runs lists them, each with its observations, its recorded inputs (t,input,value, each held from its time on), the state it started from, and where its data come from: taktwerk (a record), pilot_plant (an export), literature (digitised, with a citation), campaign (an earlier one) or synthetic. Provenance lives in the runs only; the dataset’s doc says what the runs are for:

[dataset]
contract_version = "1.0.0"
model = "refolding"
doc = "Two fed-batch refolding runs with different feeds: native and aggregated protein."

[[run]]
name = "fed-a"
source = "taktwerk"
provenance = "taktwerk record 20261019T214133Z-ef645d5c of project `fed-a` …"
observations = { path = "data/fed-a.csv", sha256 = "sha256:4b87…" }
inputs = { path = "data/fed-a.inputs.csv", sha256 = "sha256:c992…" }
initial = { v = 1.0, c_u = 0.2, c_n = 0.0, c_a = 0.0 }

Each file is pinned by its sha256; a partner’s data can stay out of git, and an estimate then asks for the file. An offline sample’s time is when it was drawn; an analysis delay belongs in the run’s doc.

Two commands add runs, dry runs until --write; the first run of a new dataset also takes --dataset-doc:

labinetix analyze import-record models/refolding/datasets/fed-refolding/record/fed-a/rows/<run-id> \
  --model-dir models/refolding --dataset fed-refolding --run fed-a --inputs f_in \
  --dataset-doc "Two fed-batch refolding runs with different feeds: native and aggregated protein." \
  --observe c_n_out=0.01,c_a_out=0.01 --every 25 \
  --project models/refolding/datasets/fed-refolding/record/fed-a.toml --write
labinetix analyze synthesize models/refolding/scenarios/fed-batch-late.toml \
  --sigma c_n_out=0.01,c_a_out=0.01 --seed 7 --at 0.333,0.715,1.172,1.618,2.084,2.547,3.013,3.468,3.927 \
  --dataset fed-refolding --run late-b --write

import-record reads a taktwerk run directory: the observed outputs every --every-th row, the recorded inputs from the tick that used them (the recorder must keep every tick), and the initial state from the project file the record pins by hash. synthesize --run records the scenario’s initial state and inputs as the run’s own.

An experiment on a dataset of runs simulates each run with the scenario’s solver and parameters, from the run’s initial state and with its recorded inputs in place of the scenario’s schedule of them, and fits all runs at once. estimate is shared by every run; estimate_per_run gives each run its own value of a parameter, estimate_initial its own initial value of a state (bounds in [initial_bounds]). A parameter that a composed model’s stoichiometric coefficient reads, such as a key yield, is refused per run: the derived coefficients follow from one value of it:

[experiment]
scenario = "fed-batch"
dataset = "fed-refolding"
estimate = ["k_f", "k_a"]
estimate_per_run = ["c_u_in"]
estimate_initial = ["c_u"]

The parameter set names the run-specific estimates c_u_in[fed-a] and c_u(0)[fed-a] in [fit], keeps the scenario’s value of a per-run parameter in [values], and adds a table per run with its observations, its share of χ² and its values:

[runs.fed-a]
n = 32
chi2 = 0.5387681782036545
parameters = { c_u_in = 10.001911356841529 }
initial = { c_u = 0.20002006164613753 }

models/refolding holds this example: two fed-batch refolding runs, one a taktwerk record of the model’s raw package, the other synthetic samples between the scenario’s steps.

Where to go from here

• models/do_probe_kla is the dynamic method of measuring kLa, with an experiment that recovers kLa, the uptake rate and the probe’s time constant, and a second scenario where a slow probe makes kLa and its lag impossible to separate.
• Composition and balances: models whose yields cannot drift out of the elemental balances, so an estimation cannot produce one that creates electrons.
• From a model to a released package: a model’s generated package under taktwerk, with a parameter set’s estimates as the instance’s values (define instance --parameter-set).
• CLI reference: every define and analyze command.