2.4. Fitting with different experimental setups
Stable-isotope labeling experiments differ in how the label is introduced and removed. The fit() and plot() methods support three experimental setups via the experiment parameter:
Mode |
Description |
Labeling curve |
|---|---|---|
|
Label on for a very long period, removed at \(t=0\) |
\(W(t) = \mathbf{s}^\top e^{\mathbf{M}t} \mathbf{1}\), decays from 1 to 0 |
|
Label on for a defined period (\(T\)), removed at \(t=0\) |
Step-up during \([-T, 0)\), step-down for \(t \ge 0\) |
|
Infinitesimally short label pulse at \(t=0\) |
\(-\dot{W}(t)\), the impulse response |
All three setups are built from the same wash-out kernel \(W(t)\), so the same compartmental model can be fitted to data from any of these designs.
This notebook demonstrates:
Generating synthetic labeling data for each experimental setup using a known two-state model.
Fitting the model to each dataset and recovering the true parameters.
Plotting the data together with the fitted labeling curves.
Visualising the fitted model as a transition-rate graph with
draw_transition_matrix().
[1]:
from symbolic_compartmental_model import SymbolicCompartmentalModel, ExperimentSetup
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path
rng = np.random.default_rng(seed=42)
# ── ground-truth parameters ──────────────────────────────────────────────────
K1_TRUE = 0.5 # self-turnover of state 0
K2_TRUE = 2.0 # self-turnover of state 1; also the contributed turnover from state 0 into state 1
A_TRUE = 0.3 # fraction of observable signal contributed by state 0 (slow pool)
NOISE = 0.03 # measurement noise standard deviation
# ── model template: cascade (state 0 → state 1) with three free parameters ──
cm = SymbolicCompartmentalModel(n_states=2)
k1 = cm.add_parameter("k1", lb=0.05, ub=5.0)
k2 = cm.add_parameter("k2", lb=0.05, ub=5.0)
a = cm.add_parameter("a", lb=0.01, ub=0.99)
cm.contributed_turnovers = [[-k1, 0.0], [k2, -k2]]
cm.observed_pool_weights = [a, 1.0 - a]
true_params = [K1_TRUE, K2_TRUE, A_TRUE]
print("Contributed-turnover matrix M:")
display(cm.contributed_turnovers)
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[1], line 1
----> 1 from symbolic_compartmental_model import SymbolicCompartmentalModel, ExperimentSetup
2 import matplotlib.pyplot as plt
3 import numpy as np
4 from pathlib import Path
ImportError: cannot import name 'ExperimentSetup' from 'symbolic_compartmental_model' (/home/docs/checkouts/readthedocs.org/user_builds/symbolic-compartmental-model/envs/stable/lib/python3.12/site-packages/symbolic_compartmental_model/__init__.py)
Each entry \(M_{ij}\) of the contributed-turnover matrix is defined as the flux between states divided by the pool size of the receiving state:
Diagonal entries \(M_{ii} \le 0\) are the self-turnover: total outflux from state \(i\) per unit pool size of state \(i\).
For the cascade matrix above:
\(M_{00} = -k_1\): self-turnover of state 0.
\(M_{10} = k_2\): the flux flowing from state 0 into state 1, divided by the pool size of state 1. This equals \(k_2\) (the self-turnover of state 1), so the row sum of row 1 is \(k_2 - k_2 = 0\): state 1 is internally mass-balanced — every molecule that clears from state 1 is replaced by one coming from state 0.
\(M_{11} = -k_2\): self-turnover of state 1.
The observed signal is a mixture of both pools: \(\mathbf{s} = (a,\, 1-a)\).
We generate synthetic data by evaluating the exact labeling predictor at a set of time points and adding Gaussian noise.
2.4.1. Washout experiment
In a washout (chase) experiment the molecule of interest is fully labeled at \(t=0\), after which labeling stops. The fraction of label remaining in the observable pool is:
This is the default mode for both fit() and plot() (experiment="washout").
[2]:
# ── generate synthetic washout data ─────────────────────────────────────────
# include early time points to capture the fast component (timescale ≈ 1/k2)
t_washout = np.array([0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 8.0])
washout_fn = cm._labeling_predictor(
x=true_params,
experiment=ExperimentSetup.WASHOUT,
period=None,
normalize=False,
use_precompiled=False,
)
y_washout = np.array([washout_fn(t) for t in t_washout])
y_washout += rng.normal(0, NOISE, size=len(y_washout))
y_washout = np.clip(y_washout, 0.0, 1.0)
# ── fit ─────────────────────────────────────────────────────────────────────
fit_washout = cm.fit(
t_washout, y_washout,
experiment=ExperimentSetup.WASHOUT,
)
print(fit_washout)
# ── plot ─────────────────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(6, 3.5), dpi=120)
ax.scatter(t_washout, y_washout, zorder=3, label="synthetic data")
fit_washout.cm.plot(
"f", ax=ax,
experiment=ExperimentSetup.WASHOUT,
label=(
f"best fit "
f"$k_1$={fit_washout.popt[0]:.2f}, "
f"$k_2$={fit_washout.popt[1]:.2f}, "
f"$a$={fit_washout.popt[2]:.2f}"
),
)
ax.axhline(0, color="black", linewidth=0.5)
ax.set_title("Washout experiment")
ax.set_xlabel("time after label removal")
ax.set_ylabel("labeling fraction $f(t)$")
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
if Path("../results").exists():
fig.savefig("../results/example_experiment_washout.svg")
display(fig)
plt.close(fig)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 3
1 # ── generate synthetic washout data ─────────────────────────────────────────
2 # include early time points to capture the fast component (timescale ≈ 1/k2)
----> 3 t_washout = np.array([0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 8.0])
4
5 washout_fn = cm._labeling_predictor(
6 x=true_params,
NameError: name 'np' is not defined
2.4.2. Step experiment
In a step experiment the label is applied for a finite period \(T\) and then removed. Time is measured from the moment labeling stops (\(t = 0\)), so measurements taken during the labeling period have negative \(t\) (as far back as \(t = -T\)):
The period argument (in the same time units as the data) tells the model how long labeling lasted.
[3]:
T_STEP = 4.0 # labeling period
# ── generate synthetic step data ──────────────────────────────────────────────
# include early post-washout points to capture the fast component
t_step = np.array([-3.0, -2.0, -1.0, 0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0])
step_fn = cm._labeling_predictor(
x=true_params,
experiment=ExperimentSetup.STEP,
period=T_STEP,
normalize=False,
use_precompiled=False,
)
y_step = np.array([step_fn(t) for t in t_step])
y_step += rng.normal(0, NOISE, size=len(y_step))
y_step = np.clip(y_step, 0.0, 1.0)
# ── fit ──────────────────────────────────────────────────────────────────────
fit_step = cm.fit(
t_step, y_step,
experiment=ExperimentSetup.STEP,
period=T_STEP,
normalize=False,
)
print(fit_step)
# ── plot ──────────────────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(6, 3.5), dpi=120)
ax.scatter(t_step, y_step, zorder=3, label="synthetic data")
fit_step.cm.plot(
"f", ax=ax,
experiment=ExperimentSetup.STEP,
period=T_STEP,
normalize=False,
label=(
f"best fit "
f"$k_1$={fit_step.popt[0]:.2f}, "
f"$k_2$={fit_step.popt[1]:.2f}, "
f"$a$={fit_step.popt[2]:.2f}"
),
)
ax.axvline(0, color="gray", linestyle="--", linewidth=1, label="label off ($t=0$)")
ax.axvline(-T_STEP, color="gray", linestyle=":", linewidth=1, label=f"label on ($t=-T$)")
ax.axhline(0, color="black", linewidth=0.5)
ax.set_title(f"Step experiment ($T = {T_STEP}$)")
ax.set_xlabel("time relative to label removal")
ax.set_ylabel("labeling fraction $f(t)$")
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
if Path("../results").exists():
fig.savefig("../results/example_experiment_step.svg")
display(fig)
plt.close(fig)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 5
1 T_STEP = 4.0 # labeling period
2
3 # ── generate synthetic step data ──────────────────────────────────────────────
4 # include early post-washout points to capture the fast component
----> 5 t_step = np.array([-3.0, -2.0, -1.0, 0.0, 0.2, 0.5, 1.0, 2.0, 3.0, 5.0])
6
7 step_fn = cm._labeling_predictor(
8 x=true_params,
NameError: name 'np' is not defined
2.4.3. Pulse experiment
In a pulse experiment an infinitesimally short bolus of label is introduced at \(t = 0\). The detected signal is the impulse response of the system — equivalently, the age distribution of molecules in the observable pool:
Because \(f(t)\) is a probability density, its values can exceed 1 — particularly for cascade models where the label must first pass through upstream compartments before reaching the observable pool. The normalization normalize=True scales the curve so that \(f(0) = 1\) rather than reporting the absolute impulse amplitude. The same flag must be set consistently in both fit() and plot().
Physical interpretation of values > 1: in the cascade below, newly synthesized molecules enter state 0 first. The observable signal grows for a while as molecules “flow through” into state 1, before eventually decaying. The peak value reflects how strongly the cascade amplifies the signal at intermediate ages.
[4]:
NOISE_PULSE = 0.04 # slightly larger noise for the pulse response
# ── generate synthetic pulse data ────────────────────────────────────────────
t_pulse = np.array([0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0])
pulse_fn = cm._labeling_predictor(
x=true_params,
experiment=ExperimentSetup.PULSE,
period=None,
normalize=True,
use_precompiled=False,
)
y_pulse = np.array([pulse_fn(t) for t in t_pulse])
y_pulse += rng.normal(0, NOISE_PULSE, size=len(y_pulse))
# normalized pulse response is a probability density: values >= 0 but can exceed 1
y_pulse = np.clip(y_pulse, 0.0, None)
# ── fit ──────────────────────────────────────────────────────────────────────
fit_pulse = cm.fit(
t_pulse, y_pulse,
experiment=ExperimentSetup.PULSE,
normalize=True,
)
print(fit_pulse)
# ── plot ──────────────────────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(6, 3.5), dpi=120)
ax.scatter(t_pulse, y_pulse, zorder=3, label="synthetic data")
fit_pulse.cm.plot(
"f", ax=ax,
experiment=ExperimentSetup.PULSE,
normalize=True,
label=(
f"best fit "
f"$k_1$={fit_pulse.popt[0]:.2f}, "
f"$k_2$={fit_pulse.popt[1]:.2f}, "
f"$a$={fit_pulse.popt[2]:.2f}"
),
)
ax.axhline(0, color="black", linewidth=0.5)
ax.set_title("Pulse experiment (normalized)")
ax.set_xlabel("time after pulse")
ax.set_ylabel("age distribution (normalized)")
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
if Path("../results").exists():
fig.savefig("../results/example_experiment_pulse.svg")
display(fig)
plt.close(fig)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 4
1 NOISE_PULSE = 0.04 # slightly larger noise for the pulse response
2
3 # ── generate synthetic pulse data ────────────────────────────────────────────
----> 4 t_pulse = np.array([0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 3.0, 5.0])
5
6 pulse_fn = cm._labeling_predictor(
7 x=true_params,
NameError: name 'np' is not defined
2.4.4. Comparing the fitted models
All three fits should recover parameters close to the ground truth (\(k_1 = 0.5\), \(k_2 = 2.0\), \(a = 0.3\)). The recovery quality depends on the number of data points, the noise level, and the identifiability of each experimental design.
We can also inspect the fitted model graphically using draw_transition_matrix(), which renders the compartment network as a directed graph. Each state node (\(s_i\)) is connected to external-input nodes (\(i_j\)) when there is net external influx, and between-state edges carry the fitted transfer rates.
[5]:
fig, axs = plt.subplots(2, 2, figsize=(8, 8), dpi=120)
for ax, params, title in zip(
axs.flat,
[true_params, fit_washout.popt, fit_step.popt, fit_pulse.popt],
["True", "Washout", "Step", "Pulse"],
):
cm.draw_transition_matrix(ax=ax, x=params)
k1_fit, k2_fit, a_fit = params
ax.set_title(
f"{title}\n $k_1$={k1_fit:.2f}, $k_2$={k2_fit:.2f}, $a$={a_fit:.2f}",
fontsize=8
)
ax.set_axis_off()
fig.tight_layout()
if Path("../results").exists():
fig.savefig("../results/example_experiment_transition_matrix.svg")
display(fig)
plt.close(fig)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 fig, axs = plt.subplots(2, 2, figsize=(8, 8), dpi=120)
2
3 for ax, params, title in zip(
4 axs.flat,
NameError: name 'plt' is not defined
Each graph shows the same two-edge cascade topology. Edge labels are contributed turnovers (flux divided by receiving pool size):
\(i_0 \to s_0\) unlabeled external input to state 0 (value = \(k_1\)), replenishing state 0’s clearance.
\(s_0 \to s_1\) contributed turnover from state 0 into state 1 (\(M_{10} = k_2\)) — the flux from state 0 into state 1, divided by the pool size of state 1.
Because state 1 is internally mass-balanced (row 1 sums to zero), there is no external-input node for state 1.
Download notebook