End-to-end system identification: a user-defined model

This notebook takes a pharmacokinetic model from a plain Python function all the way to confidence intervals on its parameters, using a user-defined model.

The point is not that the fit succeeds. It is that a good fit tells you almost nothing about whether your parameters are identifiable — and that the toolbox will tell you the difference if you ask it.

A MATLAB Live Script version (Examples/pk_user_defined.m) runs the same computation with the same section numbering, so the two can be read side by side. Throughout, # MATLAB: comments give the equivalent toolbox call.

The model

A one-compartment model with first-order absorption, the standard description of an orally administered drug:

$$c(t)=\frac{D\,k_a}{V(k_a-k_e)}\left(e^{-k_e t}-e^{-k_a t}\right)$$

Three factors are to be identified, with the dose $D=100$ mg known:

• $k_a$ — absorption rate constant (1/h)
• $k_e$ — elimination rate constant (1/h)
• $V$ — apparent volume of distribution (L)

This is the whole model. Everything below reaches it through the toolbox rather than calling it directly — model.evaluate(...), the counterpart of MATLAB's gsua_eval.

import warnings
warnings.filterwarnings("ignore")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from gsua_csb import (UserFunctionModel, design_matrix, uncertainty_analysis,
                      coverage_metric, parameter_estimation,
                      identifiability_analysis, profile_likelihood)

DOSE = 100.0        # administered dose (mg), known


def pk_absorption(params, t):
    """Plasma concentration after a single oral dose.

    params : (3,) array of [ka, ke, V]
        ka  absorption rate constant (1/h)
        ke  elimination rate constant (1/h)
        V   apparent volume of distribution (L)
    t : (n,) array of sampling times, in hours

    Returns
    -------
    (1, n) array -- one observed output over the n times. The leading axis is the
    output index: gsua_csb models always return (n_outputs, n_times), even for a
    single output, which is what makes `model.output` selection possible.

    MATLAB equivalent: Examples/pkAbsorptionModel.m.
    """
    ka, ke, V = params
    return ((DOSE * ka) / (V * (ka - ke)) * (np.exp(-ke * t) - np.exp(-ka * t)))[None, :]

1. Preparing the environment

UserFunctionModel wraps the function above into the object every other gsua_csb routine consumes — the counterpart of the summary table T that gsua_dataprep builds in MATLAB.

The $k_a$ bounds are deliberately kept above the $k_e$ bounds. At $k_a=k_e$ the closed form is singular, and swapping the two leaves $c(t)$ unchanged — the classic flip-flop ambiguity. Excluding it keeps this example about experimental design rather than an algebraic accident.

# Factor bounds: one row per factor, [lower, upper] -- same layout as MATLAB's `ranges`.
ranges = np.array([[0.6,  3.0],     # ka  absorption rate (1/h)
                   [0.05, 0.5],     # ke  elimination rate (1/h)
                   [5.0,  40.0]])   # V   volume of distribution (L)

truth = np.array([1.2, 0.25, 15.0])                        # values to recover
xdata = np.array([0.25, 0.5, 1, 1.5, 2, 3, 4, 6, 8,
                  10, 12, 16, 20, 24.0])                   # sampling schedule (hours)

# MATLAB: [T,~] = gsua_dataprep('pkAbsorptionModel', ranges, 'domain',[0 24], ...)
model = UserFunctionModel(
    func=pk_absorption,          # <-- the model function defined above
    names=["ka", "ke", "V"],
    range=ranges,
    nominal=truth,
    domain=xdata,
    output_names=["concentration"],
)
print("factors:", model.names, "| free:", int((~model.fixed).sum()))
print("outputs:", model.active_output_names)

factors: ['ka', 'ke', 'V'] | free: 3
outputs: ['concentration']

2. Synthetic data

Working from synthetic data means the truth is known, so the confidence intervals can be checked rather than merely reported.

Note that the model is evaluated through model.evaluate, not by calling pk_absorption directly. That is the toolbox's evaluation path — it applies the model's output selection (see model.output) and matches what every routine below does internally. MATLAB's gsua_eval plays the same role.

rng = np.random.default_rng(0)                  # fix the noise draw so the page reproduces

# MATLAB: clean = gsua_eval(truth, T, xdata, [], false, false);
clean = model.evaluate(truth, xdata)            # noise-free model output at the sample times
ydata = clean * (1 + 0.08 * rng.standard_normal(clean.shape))   # 8% proportional noise

tdense = np.linspace(0.05, 24, 300)             # dense grid, for drawing curves only

plt.figure(figsize=(6.4, 3.8))
plt.plot(tdense, model.evaluate(truth, tdense)[0], lw=1.6, label="true model")
plt.plot(xdata, ydata[0], "ko", ms=5, label="measurements")
plt.xlabel("time (h)"); plt.ylabel("concentration (mg/L)")
plt.title("Simulated single-dose concentration data")
plt.legend(); plt.grid(alpha=.3); plt.tight_layout(); plt.show()

3. Can these parameters be estimated at all?

Before spending any optimizer budget it is worth asking whether the data is even reachable: does it fall inside the range of behaviours the model can produce over the bounds declared in section 1? If it does not, no amount of optimization will help — the model structure or the bounds are wrong, and that has to be fixed first. This is the reachability check that opens the toolbox's semi-automated identification cycle.

design_matrix samples the factor box and uncertainty_analysis runs the Monte-Carlo ensemble over those samples.

# MATLAB: M0 = gsua_dmatrix(T,500);  Y0 = gsua_ua(M0,T,'xdata',xdata,'ynom',ydata);
M0 = design_matrix(model, 500, seed=0)          # 500 samples of the factor box
ua0 = uncertainty_analysis(model, M0, xdata=xdata, y_exp=ydata)

# coverage_metric reduces the ensemble to a 5-95% band and scores it.
# MATLAB: [cost_data,cost_band,P5,~,P95] = gsua_covmetric(Y0, ydata, 'margin',0.1);
cost_data, cost_band, p5, _, p95 = coverage_metric(ua0.Y, ydata, margin=0.1)
contained = float(np.mean((ydata >= p5) & (ydata <= p95)))

pd.DataFrame({"contained": [contained],
              "median_band_width": [float(np.median(p95 - p5))],
              "cost_data": [cost_data],
              "cost_band": [cost_band]}).round(4)

	contained	median_band_width	cost_data	cost_band
0	1.0	4.8627	127.5652	2.246380e+07

All of the data lies inside the reachable band, so estimation is worth attempting.

cost_data and cost_band are both far above 1 here, and that is expected rather than alarming: they are normalized against a tight tolerance and are meaningful as a post-convergence check (section 7), not as a pass/fail gate against a prior range this wide. The containment fraction is the number to read at this stage. The band is enormous — a median width of several mg/L against data that never exceeds 4.5 — which is what an uninformative prior looks like before any fitting.

plt.figure(figsize=(6.4, 3.8))
plt.plot(xdata, p5[0], color="0.4", lw=1, label="5th percentile")
plt.plot(xdata, p95[0], color="0.4", lw=1, label="95th percentile")
plt.fill_between(xdata, p5[0], p95[0], color="0.85")
plt.plot(xdata, ydata[0], "ko", ms=5, label="measurements")
plt.xlabel("time (h)"); plt.ylabel("concentration (mg/L)")
plt.title(f"Reachable band before fitting ({contained:.0%} of data contained)")
plt.legend(); plt.grid(alpha=.3); plt.tight_layout(); plt.show()

4. Estimating the parameters

parameter_estimation runs a multistart fit: n=20 restarts the optimizer from twenty different points in the factor space, which is how you find out whether the problem has one solution or several.

# margin=0.1 selects the correlation-penalized cost and records the margin on the
# result, so functions further down the pipeline can recover what was scored.
# MATLAB: [T3,res3] = gsua_pe(T, xdata, ydata, 'solver','lsqc', 'N',20, 'margin',0.1)
pe3 = parameter_estimation(model, xdata, ydata, n=20,
                           solver="least_squares", margin=0.1, seed=0)

best3 = pe3.x[np.argmin(pe3.cost)]              # pe3.x is (20, 3): one row per restart
pd.DataFrame({"true": truth, "estimated": best3.round(4)}, index=model.names)

	true	estimated
ka	1.20	1.1554
ke	0.25	0.2507
V	15.00	14.5982

print(f"cost across the 20 restarts:  min {pe3.cost.min():.5g}   max {pe3.cost.max():.5g}")

plt.figure(figsize=(6.4, 3.8))
plt.plot(tdense, model.evaluate(best3, tdense)[0], lw=1.6, label="fitted model")
plt.plot(xdata, ydata[0], "ko", ms=5, label="measurements")
plt.xlabel("time (h)"); plt.ylabel("concentration (mg/L)")
plt.title(f"Fit with all three factors free (cost = {pe3.cost.min():.4g})")
plt.legend(); plt.grid(alpha=.3); plt.tight_layout(); plt.show()

cost across the 20 restarts:  min 0.18284   max 0.18284

Every one of the twenty restarts converged to the same cost, and the fitted curve passes cleanly through the data. On most projects this is where the analysis would stop.

5. Diagnosing identifiability

The correlation between the repeated estimates is the first warning sign. Values near $\pm 1$ mean the factors trade off against each other: many different combinations reproduce the same curve.

# MATLAB: array2table(corr(Est3'), ...)
ia3 = identifiability_analysis(model, pe3.x, cost=pe3.cost, cost_rtol=0.1, seed=0)
pd.DataFrame(ia3.correlation.round(4), index=model.names, columns=model.names)

	ka	ke	V
ka	1.0000	-0.9999	1.0
ke	-0.9999	1.0000	-1.0
V	1.0000	-1.0000	1.0

# profile_likelihood steps each factor away from the estimate while re-fitting all the
# others at every step, and reports where the fit stays statistically acceptable.
# margin is the assumed relative standard deviation: 0.08 matches the 8% noise above.
# MATLAB's gsua_likelihood offsets this by one, so the same run there is margin=1.08:
# MATLAB: gsua_likelihood(T3, xdata, ydata, 0.95, 0.05, 1.08, 0.01, 0.01, 15, 1, ...)
pl3 = profile_likelihood(model, xdata, ydata, alpha=0.95, margin=0.08)

pd.DataFrame({
    "CI_low":     pl3.range[:, 0].round(4),
    "CI_high":    pl3.range[:, 1].round(4),
    "width":      (pl3.range[:, 1] - pl3.range[:, 0]).round(4),
    "prior_low":  ranges[:, 0],
    "prior_high": ranges[:, 1],
}, index=model.names)

	CI_low	CI_high	width	prior_low	prior_high
ka	0.9352	1.4798	0.5447	0.60	3.0
ke	0.2448	0.2709	0.0261	0.05	0.5
V	13.0391	16.4258	3.3867	5.00	40.0

This is the result worth stopping on. The correlation between $k_a$ and $k_e$ is essentially $-1$: the three factors cannot be separated from a single oral concentration curve, however well that curve fits. That is a textbook pharmacokinetic result, not a failure of the optimizer.

6. The remedy: fix what another experiment already knows

The standard resolution is to measure $V$ separately, in an intravenous study where it is directly identifiable, and then estimate only the two rate constants. A factor is fixed by giving it a degenerate range — lower bound equal to upper bound — which is what Model.fix does in place.

import copy
model_fixed = copy.deepcopy(model)
model_fixed.fix("V", 15.0)                      # MATLAB: give the V row the range [15 15]
free = np.where(~model_fixed.fixed)[0]          # -> [0, 1], i.e. ka and ke
print("free factors now:", [model.names[i] for i in free])

free factors now: ['ka', 'ke']

pe2 = parameter_estimation(model_fixed, xdata, ydata, n=20,
                           solver="least_squares", margin=0.1, seed=0)
best2 = pe2.x[np.argmin(pe2.cost)]
ia2 = identifiability_analysis(model_fixed, pe2.x, cost=pe2.cost, cost_rtol=0.1, seed=0)
pl2 = profile_likelihood(model_fixed, xdata, ydata, alpha=0.95, margin=0.08,
                         params=list(free))    # profile only the free factors

pd.DataFrame({
    "true":      truth[free],
    "estimated": best2[free].round(4),
    "CI_low":    pl2.range[free, 0].round(4),
    "CI_high":   pl2.range[free, 1].round(4),
    "width":     (pl2.range[free, 1] - pl2.range[free, 0]).round(4),
}, index=[model.names[i] for i in free])

	true	estimated	CI_low	CI_high	width
ka	1.20	1.2082	1.0064	1.4433	0.4369
ke	0.25	0.2418	0.2480	0.2619	0.0139

Both remaining factors are now recovered close to the truth, and both intervals sit well inside their bounds instead of running to them.

7. The two runs side by side

pd.DataFrame({
    "all three free": [pe3.cost.min(), ia3.correlation[0, 1],
                       pl3.range[0, 1] - pl3.range[0, 0],
                       pl3.range[1, 1] - pl3.range[1, 0]],
    "V fixed":        [pe2.cost.min(), ia2.correlation[0, 1],
                       pl2.range[0, 1] - pl2.range[0, 0],
                       pl2.range[1, 1] - pl2.range[1, 0]],
}, index=["best cost", "corr(ka, ke)", "CI width ka", "CI width ke"]).round(4)

	all three free	V fixed
best cost	0.1828	0.1899
corr(ka, ke)	-0.9999	0.3869
CI width ka	0.5447	0.4369
CI width ke	0.0261	0.0139

x = np.arange(2); w = 0.35
plt.figure(figsize=(5.6, 3.6))
plt.bar(x - w/2, [pl3.range[i, 1] - pl3.range[i, 0] for i in range(2)], w, label="all three free")
plt.bar(x + w/2, [pl2.range[i, 1] - pl2.range[i, 0] for i in range(2)], w, label="V fixed")
plt.xticks(x, ["ka", "ke"]); plt.ylabel("95% confidence interval width")
plt.title("Fixing one factor sharpens the other two")
plt.legend(); plt.grid(alpha=.3, axis="y"); plt.tight_layout(); plt.show()

8. What this example shows

Fixing $V$ made the fit slightly worse and the science considerably better: the correlation between $k_a$ and $k_e$ collapses, both intervals tighten, and the estimates move onto the truth.

Cost measures how well a curve passes through points. It does not measure whether the factors that produced that curve could have been recovered. Only the identifiability analysis answers that, which is why it belongs inside the workflow rather than after it.

Comparing with the MATLAB version. Both implementations use the same acceptance threshold (chi-squared at one degree of freedom, halved) — but margin is offset by one between them. MATLAB's gsua_likelihood treats margin-1 as the assumed relative standard deviation, so its margin=1.08 is this notebook's margin=0.08; both mean 8%. Passing 0.1 to MATLAB would assert 90% noise, and would also flip gsua_pe's internal +1 offset positive, silently switching its inner refit from the likelihood to plain least squares. With the conventions matched the two agree closely — $k_a$ interval widths of 0.48 (MATLAB) and 0.61 (Python) on this problem.

The companion symbolic-ODE example reaches the same conclusion from the opposite direction: there, the dataset that fits better is the one whose parameters are less identifiable.

Where a multistart run does spread across the factor space instead of converging to a single point, identifiability_analysis adds detection of multiple global minima, noise_floor calibrates which fits to accept against the observation noise, and design_matrix(..., method="joint") propagates the accepted set without destroying its correlation structure.