43.4 Building a Calibrated Simulator

"I calibrated my model, then asked it what would happen if I doubled the degradation rate. It told me. Then I checked in the lab. It was right. That was the moment I stopped calling it a toy."

A Posterior Distribution That Earned Its Keep

1. The Biological System: A Genetic Toggle Switch

Imagine a cell carrying just a few dozen molecules of two rival proteins, each trying to silence the other; in the next hour, random chance will tip the balance and lock the cell into one of two fates, yet no one can predict which fate by watching averages alone.

The genetic toggle switch, first constructed synthetically by Gardner et al. (2000), consists of two genes that mutually repress each other. Gene A produces protein A, which represses gene B. Gene B produces protein B, which represses gene A. This mutual inhibition creates two stable steady states: high A / low B, or low A / high B. Stochastic fluctuations in molecule counts can cause the system to switch between states, a phenomenon with implications for cell fate decisions, antibiotic resistance, and developmental biology.

The reaction network has six reactions:

The Hill coefficient \(n\) controls the cooperativity of repression (the degree to which multiple repressor molecules must bind together to shut down a gene, producing an all-or-nothing response rather than a gradual one): higher \(n\) produces sharper switches and more stable bistability. The production rates \(\alpha_A, \alpha_B\) and degradation rates \(\delta_A, \delta_B\) are the parameters to calibrate. In short: a calibrated simulator turns "I think this drug will work" into "the posterior probability that this drug breaks the switch is 0.83," replacing intuition with quantified confidence.

The full calibration pipeline, illustrated in Figure 43.4.1, connects these components into a four-stage workflow: simulate, compare, refine, and interrogate.

2. Implementing the Toggle Switch Simulator

import numpy as np
from dataclasses import dataclass


@dataclass
class ToggleSwitchParams:
    """Parameters for the genetic toggle switch model."""
    alpha_a: float = 50.0   # max production rate of A
    alpha_b: float = 50.0   # max production rate of B
    delta_a: float = 1.0    # degradation rate of A
    delta_b: float = 1.0    # degradation rate of B
    n: float = 2.0          # Hill coefficient (cooperativity)


def toggle_switch_gillespie(params, t_max=100.0, seed=42):
    """Simulate the genetic toggle switch using the Gillespie algorithm.

    State: [A, B] (molecule counts of the two proteins)

    Args:
        params: ToggleSwitchParams or dict with parameter values
        t_max: simulation duration
        seed: random seed

    Returns:
        dict with 'times' and 'states' arrays
    """
    if isinstance(params, dict):
        p = ToggleSwitchParams(**params)
    else:
        p = params

    rng = np.random.default_rng(seed)

    # Initial state: start near one steady state
    x = np.array([0.0, 0.0])  # [A, B]
    t = 0.0

    times = [t]
    states = [x.copy()]

    while t < t_max:
        A, B = x[0], x[1]

        # Propensities with Hill-function repression
        a1 = p.alpha_a / (1.0 + B**p.n)    # production of A
        a2 = p.delta_a * A                  # degradation of A
        a3 = p.alpha_b / (1.0 + A**p.n)    # production of B
        a4 = p.delta_b * B                  # degradation of B

        propensities = np.array([a1, a2, a3, a4])
        a0 = propensities.sum()

        if a0 == 0:
            break

        # Time to next event
        tau = rng.exponential(1.0 / a0)
        t += tau
        if t > t_max:
            break

        # Which reaction fires?
        j = rng.choice(4, p=propensities / a0)

        # Stoichiometry: [dA, dB]
        stoich = np.array([
            [1, 0],    # R1: produce A
            [-1, 0],   # R2: degrade A
            [0, 1],    # R3: produce B
            [0, -1],   # R4: degrade B
        ])

        x += stoich[j]
        x = np.maximum(x, 0)

        times.append(t)
        states.append(x.copy())

    return {"times": np.array(times), "states": np.array(states)}


# Run a single trajectory
params = ToggleSwitchParams(alpha_a=50, alpha_b=50, delta_a=1, delta_b=1, n=2)
result = toggle_switch_gillespie(params, t_max=200.0, seed=42)
print(f"Toggle switch: {len(result['times'])} events over "
      f"{result['times'][-1]:.1f} time units")
print(f"Final state: A={result['states'][-1, 0]:.0f}, "
      f"B={result['states'][-1, 1]:.0f}")

3. Generating Synthetic Experimental Data

Real discovery workflows use experimental measurements such as single-cell fluorescence time series. Here, synthetic data generated from known ground-truth parameters stands in for real observations, allowing verification that calibration recovers the true values before the pipeline meets real data.

def generate_observed_data(
    true_params, n_cells=20, t_max=200.0,
    measurement_times=None, noise_std=5.0, seed=42,
):
    """Generate synthetic experimental data from the toggle switch model.

    Simulates n_cells independent cells, samples at measurement_times,
    and adds Gaussian measurement noise.

    Args:
        true_params: ToggleSwitchParams with ground-truth values
        n_cells: number of independent cell trajectories
        t_max: simulation duration
        measurement_times: array of time points to sample
        noise_std: standard deviation of measurement noise
        seed: random seed

    Returns:
        observed: dict with measurement_times, cell_data, and summaries
    """
    rng = np.random.default_rng(seed)

    if measurement_times is None:
        measurement_times = np.linspace(10, t_max, 20)

    all_measurements = []

    for i in range(n_cells):
        cell_seed = rng.integers(0, 2**31)
        result = toggle_switch_gillespie(true_params, t_max=t_max, seed=cell_seed)

        # Interpolate to measurement times
        A_measured = np.interp(
            measurement_times, result["times"], result["states"][:, 0]
        )
        B_measured = np.interp(
            measurement_times, result["times"], result["states"][:, 1]
        )

        # Add measurement noise
        A_noisy = A_measured + rng.normal(0, noise_std, len(measurement_times))
        B_noisy = B_measured + rng.normal(0, noise_std, len(measurement_times))

        all_measurements.append(np.column_stack([A_noisy, B_noisy]))

    cell_data = np.array(all_measurements)  # (n_cells, n_times, 2)

    return {
        "measurement_times": measurement_times,
        "cell_data": cell_data,
        "n_cells": n_cells,
    }


# Ground truth parameters
true_params = ToggleSwitchParams(
    alpha_a=40.0, alpha_b=45.0, delta_a=0.8, delta_b=1.2, n=2.5
)

observed = generate_observed_data(true_params, n_cells=20, seed=42)
print(f"Generated data: {observed['n_cells']} cells, "
      f"{len(observed['measurement_times'])} time points each")

Synthetic data gives us a target to calibrate against, but comparing raw time series between simulator output and observations is impractical when trajectories are stochastic and high-dimensional; we need compact numerical summaries that capture the system's essential behavior.

4. Designing Summary Statistics for the Toggle Switch

The toggle switch has distinctive statistical signatures that good summary statistics must capture: bimodal distributions (where the histogram of protein levels shows two distinct peaks, one for each stable state), anti-correlation between A and B, and specific steady-state levels.

def toggle_switch_summaries(sim_result, measurement_times=None):
    """Compute summary statistics for toggle switch calibration.

    Extracts features that capture the bistable dynamics:
    - Mean and variance of each protein across cells/time
    - Bimodality coefficient (distinguishes unimodal from bimodal)
    - Cross-correlation between A and B
    - Fraction of time spent in each state

    Args:
        sim_result: dict with 'times' and 'states' from the simulator
        measurement_times: optional array of time points to evaluate

    Returns:
        1D numpy array of summary statistics
    """
    times = sim_result["times"]
    states = sim_result["states"]
    A = states[:, 0]
    B = states[:, 1]

    # Use the second half of the trajectory (after transient)
    midpoint = len(times) // 2
    A_ss = A[midpoint:]
    B_ss = B[midpoint:]

    summaries = []

    # Moments for each protein
    for series, name in [(A_ss, "A"), (B_ss, "B")]:
        summaries.append(np.mean(series))
        summaries.append(np.std(series))
        summaries.append(np.median(series))

        # Bimodality coefficient: (skewness^2 + 1) / kurtosis
        # Values > 5/9 suggest bimodality
        if np.std(series) > 0:
            from scipy.stats import skew, kurtosis
            s = skew(series)
            k = kurtosis(series, fisher=False)  # excess=False -> Pearson
            bimod = (s**2 + 1) / k if k > 0 else 0
            summaries.append(bimod)
        else:
            summaries.append(0.0)

    # Cross-correlation
    if np.std(A_ss) > 0 and np.std(B_ss) > 0:
        min_len = min(len(A_ss), len(B_ss))
        corr = np.corrcoef(A_ss[:min_len], B_ss[:min_len])[0, 1]
        summaries.append(corr if np.isfinite(corr) else 0.0)
    else:
        summaries.append(0.0)

    # Fraction of time in high-A state (A > B)
    min_len = min(len(A_ss), len(B_ss))
    frac_high_a = np.mean(A_ss[:min_len] > B_ss[:min_len])
    summaries.append(frac_high_a)

    # Mean switching frequency (state changes per unit time)
    if min_len > 1:
        state_labels = (A_ss[:min_len] > B_ss[:min_len]).astype(int)
        switches = np.sum(np.abs(np.diff(state_labels)))
        dt = times[-1] - times[midpoint] if times[-1] > times[midpoint] else 1.0
        summaries.append(switches / dt)
    else:
        summaries.append(0.0)

    return np.array(summaries)


def observed_data_summaries(observed):
    """Compute summary statistics from the observed multi-cell data.

    Aggregates across cells to produce population-level summaries
    that match the single-trajectory summaries from the simulator.
    """
    cell_data = observed["cell_data"]  # (n_cells, n_times, 2)

    # Average across cells, then compute summaries
    mean_trajectory = cell_data.mean(axis=0)  # (n_times, 2)

    # Use second half
    mid = len(mean_trajectory) // 2
    A_ss = mean_trajectory[mid:, 0]
    B_ss = mean_trajectory[mid:, 1]

    summaries = []

    for series in [A_ss, B_ss]:
        summaries.append(np.mean(series))
        summaries.append(np.std(series))
        summaries.append(np.median(series))
        # Bimodality from cell-level distribution at the last time point
        if np.std(series) > 0:
            from scipy.stats import skew, kurtosis
            s = skew(series)
            k = kurtosis(series, fisher=False)
            summaries.append((s**2 + 1) / k if k > 0 else 0)
        else:
            summaries.append(0.0)

    if np.std(A_ss) > 0 and np.std(B_ss) > 0:
        summaries.append(np.corrcoef(A_ss, B_ss)[0, 1])
    else:
        summaries.append(0.0)

    summaries.append(np.mean(A_ss > B_ss))

    # Switching frequency from individual cells
    switch_rates = []
    for cell in cell_data:
        mid_c = len(cell) // 2
        labels = (cell[mid_c:, 0] > cell[mid_c:, 1]).astype(int)
        switch_rates.append(np.sum(np.abs(np.diff(labels))))
    summaries.append(np.mean(switch_rates))

    return np.array(summaries)


obs_summaries = observed_data_summaries(observed)
print(f"Observed summary statistics ({len(obs_summaries)} dimensions):")
print(obs_summaries.round(3))

5. Running ABC-SMC Calibration

The simulator, observed data, and summary statistics are now in place for ABC-SMC calibration. We use pyABC for its adaptive epsilon scheduling and efficient implementation.

Without principled calibration, a simulator's predictions are little more than guesswork dressed in code. Teams that skip this step routinely discover, after months of wet-lab follow-up, that their in silico drug screen was driven by a plausible but wrong parameter regime, wasting reagents and time on interventions the model never truly endorsed.

ABC-SMC is a likelihood-free inference algorithm that estimates the posterior distribution over model parameters when the likelihood function is too complex to evaluate directly, as is the case for stochastic simulators. It matters because most interesting simulators (including Gillespie models) produce stochastic outputs whose probability cannot be written in closed form, making standard Bayesian methods like Markov chain Monte Carlo (MCMC) inapplicable. The algorithm proposes candidate parameters from a prior, simulates data from each candidate, and compares the simulated summary statistics to the observed ones via a distance threshold (epsilon). It then tightens that threshold across successive populations so that surviving particles (where each particle is a candidate parameter vector that has been accepted by the distance criterion) converge toward the true posterior. The convergence relies on importance weighting: at each new population, particles from the previous round are perturbed by a kernel (typically a Gaussian centered on the old particle), and weights correct for the mismatch between this proposal distribution and the prior, ensuring the final sample approximates the true posterior rather than the proposal. Use ABC-SMC when your simulator is a black box with no tractable likelihood; if your model does have a closed-form likelihood, standard MCMC (e.g., Hamiltonian Monte Carlo via Stan or PyMC; as of 2024, PyMC v5 and NumPyro are the most actively maintained options) will be more sample-efficient.

import pyabc
from pyabc import ABCSMC, RV, Distribution, MedianEpsilon


def run_toggle_switch_calibration(
    observed, n_populations=8, population_size=300, seed=42,
):
    """Calibrate the toggle switch model using pyABC's ABC-SMC.

    Args:
        observed: dict from generate_observed_data
        n_populations: number of SMC generations
        population_size: particles per generation
        seed: random seed

    Returns:
        history: pyABC History object
        posterior_samples: DataFrame of accepted parameters
    """
    obs_summ = observed_data_summaries(observed)

    # Prior distributions (broad, covering plausible biological ranges)
    prior = Distribution(
        alpha_a=RV("uniform", 10, 90),    # 10 to 100
        alpha_b=RV("uniform", 10, 90),    # 10 to 100
        delta_a=RV("uniform", 0.1, 2.9),  # 0.1 to 3.0
        delta_b=RV("uniform", 0.1, 2.9),  # 0.1 to 3.0
        n=RV("uniform", 1.0, 4.0),        # 1.0 to 5.0
    )

    def model(params_dict):
        """Run the simulator and return summary statistics as a dict."""
        p = ToggleSwitchParams(**params_dict)
        result = toggle_switch_gillespie(p, t_max=200.0, seed=seed)
        summaries = toggle_switch_summaries(result)
        return {f"s{i}": float(v) for i, v in enumerate(summaries)}

    # Convert observed summaries to matching dict format
    obs_dict = {f"s{i}": float(v) for i, v in enumerate(obs_summ)}

    # Distance and epsilon
    distance = pyabc.PNormDistance(p=2)
    epsilon = MedianEpsilon(initial_epsilon=100.0)

    abc = ABCSMC(
        models=model,
        parameter_priors=prior,
        distance_function=distance,
        population_size=population_size,
        eps=epsilon,
    )

    import tempfile, os
    db_path = os.path.join(tempfile.gettempdir(), "toggle_calibration.db")
    abc.new(f"sqlite:///{db_path}", obs_dict)

    history = abc.run(
        minimum_epsilon=1.0,
        max_nr_populations=n_populations,
    )

    # Extract posterior samples from the final generation
    df, w = history.get_distribution()
    print("\nPosterior summary:")
    print(df.describe().round(3))
    print(f"\nTrue values: alpha_a=40, alpha_b=45, delta_a=0.8, "
          f"delta_b=1.2, n=2.5")

    return history, df


# Run calibration (this takes several minutes)
# history, posterior = run_toggle_switch_calibration(observed)

6. Validating the Calibration

Before using a calibrated simulator for counterfactual predictions, we must verify that the inferred posterior actually produces simulations consistent with the observed data. This step, called a posterior predictive check (a validation procedure that simulates new data from the inferred parameter posterior and compares it to the original observations, testing whether the calibrated model can reproduce the data it was trained on), closes a gap that calibration alone leaves open: ABC-SMC finds parameters whose summary statistics fall within a distance threshold of the observations, but it does not guarantee that the full simulated trajectories look realistic. A posterior predictive check evaluates that broader consistency. Section 7 implements this check as the validate() method of the pipeline class.

7. Counterfactual Interventions

A calibrated simulator is not just a fitted model; it is an experimental platform. You can ask "What if?" questions that are impossible, expensive, or unethical to answer in the lab. These counterfactual interventions (experiments that modify one or more model parameters to simulate a hypothetical scenario that did not actually occur) are the primary scientific payoff of calibration.

For the toggle switch, mechanistically interesting counterfactuals include:

• What if we doubled the degradation rate of protein A? (Mimics adding a protease.)
• What if we reduced the Hill coefficient? (Mimics weakening cooperative binding.)
• What if we broke the symmetry by increasing production of B? (Mimics an environmental signal.)

def counterfactual_experiment(
    posterior_samples, intervention_fn, n_posterior_draws=50,
    n_replicates=10, t_max=200.0, seed=42,
):
    """Run a counterfactual experiment using the calibrated posterior.

    For each posterior sample, applies the intervention and simulates,
    then compares to the unperturbed simulation.

    Args:
        posterior_samples: DataFrame with columns for each parameter
        intervention_fn: callable(params_dict) -> modified params_dict
        n_posterior_draws: number of posterior samples to use
        n_replicates: simulations per posterior sample
        t_max: simulation duration
        seed: random seed

    Returns:
        baseline_summaries: summary stats without intervention
        intervention_summaries: summary stats with intervention
    """
    rng = np.random.default_rng(seed)

    # Sample from the posterior
    if len(posterior_samples) > n_posterior_draws:
        indices = rng.choice(
            len(posterior_samples), n_posterior_draws, replace=False
        )
        samples = posterior_samples.iloc[indices]
    else:
        samples = posterior_samples

    baseline_results = []
    intervention_results = []

    for _, row in samples.iterrows():
        params_dict = row.to_dict()

        for rep in range(n_replicates):
            run_seed = rng.integers(0, 2**31)

            # Baseline (no intervention)
            p_base = ToggleSwitchParams(**params_dict)
            result_base = toggle_switch_gillespie(
                p_base, t_max=t_max, seed=run_seed
            )
            baseline_results.append(
                toggle_switch_summaries(result_base)
            )

            # Intervention
            modified_dict = intervention_fn(params_dict.copy())
            p_mod = ToggleSwitchParams(**modified_dict)
            result_mod = toggle_switch_gillespie(
                p_mod, t_max=t_max, seed=run_seed
            )
            intervention_results.append(
                toggle_switch_summaries(result_mod)
            )

    return np.array(baseline_results), np.array(intervention_results)


def analyze_counterfactual(baseline, intervention, stat_names=None):
    """Compare baseline and intervention summary statistics.

    Computes effect sizes and credible intervals for each statistic.
    """
    if stat_names is None:
        stat_names = [f"stat_{i}" for i in range(baseline.shape[1])]

    print(f"{'Statistic':<25} {'Baseline':>10} {'Intervention':>12} "
          f"{'Effect':>10} {'95% CI':>20}")
    print("-" * 80)

    for j, name in enumerate(stat_names):
        base_mean = np.mean(baseline[:, j])
        int_mean = np.mean(intervention[:, j])
        effect = int_mean - base_mean

        # Bootstrap 95% CI for the effect
        effects = intervention[:, j] - baseline[:, j]
        ci_lo = np.percentile(effects, 2.5)
        ci_hi = np.percentile(effects, 97.5)

        print(f"{name:<25} {base_mean:>10.2f} {int_mean:>12.2f} "
              f"{effect:>+10.2f} [{ci_lo:>+8.2f}, {ci_hi:>+8.2f}]")


# Example: What happens if we double the degradation of protein A?
def double_delta_a(params):
    params["delta_a"] *= 2.0
    return params

# Example: What happens if we reduce the Hill coefficient?
def reduce_cooperativity(params):
    params["n"] = max(1.0, params["n"] - 1.0)
    return params

stat_names = [
    "mean_A", "std_A", "median_A", "bimod_A",
    "mean_B", "std_B", "median_B", "bimod_B",
    "corr_AB", "frac_high_A", "switch_freq",
]

# In practice, after calibration:
# baseline, intervened = counterfactual_experiment(
#     posterior, double_delta_a, n_posterior_draws=50, n_replicates=10
# )
# analyze_counterfactual(baseline, intervened, stat_names)

With calibration validated and counterfactual experiments demonstrating the model's predictive power, the remaining step is to package these components into a single reusable pipeline that other researchers (and automated systems) can invoke without reassembling the pieces by hand.

8. The Complete Pipeline

The entire workflow, following the stages shown in Figure 43.4.1, assembles into a single, reusable pipeline class that integrates with the Discovery Workbench architecture from Chapter 6.

from dataclasses import dataclass, field
from typing import Callable, Optional
import json


@dataclass
class CalibratedSimulatorPipeline:
    """End-to-end pipeline for simulation-based discovery.

    Stages:
    1. Define the simulator and summary statistics
    2. Calibrate against observed data using ABC-SMC
    3. Validate the calibration (posterior predictive check)
    4. Run counterfactual experiments
    5. Report results with full provenance
    """

    simulator_fn: Callable
    summary_fn: Callable
    param_names: list
    prior_ranges: dict
    observed_data: Optional[dict] = None
    posterior_samples: Optional[np.ndarray] = None
    calibration_history: Optional[object] = None
    provenance: dict = field(default_factory=dict)

    def set_observed_data(self, data):
        """Register the observed experimental data."""
        self.observed_data = data
        self.provenance["data_registered"] = True
        self.provenance["n_observations"] = (
            data.get("n_cells", "unknown")
        )

    def calibrate(self, n_populations=8, population_size=300, seed=42):
        """Run ABC-SMC calibration."""
        if self.observed_data is None:
            raise ValueError("Call set_observed_data() first")

        obs_summaries = self.summary_fn(self.observed_data)

        # ABC-SMC calibration (using the manual implementation for clarity)
        from scipy.stats import uniform

        def prior_sampler(rng):
            params = []
            for name in self.param_names:
                lo, hi = self.prior_ranges[name]
                params.append(rng.uniform(lo, hi))
            return np.array(params)

        def prior_pdf(theta):
            for j, name in enumerate(self.param_names):
                lo, hi = self.prior_ranges[name]
                if theta[j] < lo or theta[j] > hi:
                    return 0.0
            pdf = 1.0
            for j, name in enumerate(self.param_names):
                lo, hi = self.prior_ranges[name]
                pdf *= 1.0 / (hi - lo)
            return pdf

        def sim_wrapper(theta, seed=42):
            params_dict = dict(zip(self.param_names, theta))
            return self.simulator_fn(params_dict, seed=seed)

        def summary_wrapper(sim_data):
            return self.summary_fn(sim_data)

        # Use decreasing epsilon schedule
        epsilons = np.geomspace(100, 2, n_populations)

        from scipy.spatial.distance import euclidean

        particles, weights, populations = abc_smc(
            simulator_fn=sim_wrapper,
            summary_fn=summary_wrapper,
            observed_summaries=obs_summaries,
            prior_sampler=prior_sampler,
            prior_pdf=prior_pdf,
            epsilons=epsilons,
            n_particles=population_size,
            seed=seed,
        )

        self.posterior_samples = particles
        self.provenance["calibration"] = {
            "method": "ABC-SMC",
            "n_populations": n_populations,
            "population_size": population_size,
            "final_epsilon": float(epsilons[-1]),
            "seed": seed,
        }

        return particles, weights

    def validate(self, n_checks=100, seed=42):
        """Run posterior predictive check."""
        if self.posterior_samples is None:
            raise ValueError("Call calibrate() first")

        obs_summaries = self.summary_fn(self.observed_data)
        rng = np.random.default_rng(seed)
        distances = []

        indices = rng.choice(
            len(self.posterior_samples), min(n_checks, len(self.posterior_samples)),
            replace=False,
        )

        for idx in indices:
            theta = self.posterior_samples[idx]
            params_dict = dict(zip(self.param_names, theta))
            sim_data = self.simulator_fn(
                params_dict, seed=rng.integers(0, 2**31)
            )
            sim_summ = self.summary_fn(sim_data)
            from scipy.spatial.distance import euclidean
            distances.append(euclidean(sim_summ, obs_summaries))

        distances = np.array(distances)
        self.provenance["validation"] = {
            "median_distance": float(np.median(distances)),
            "p95_distance": float(np.percentile(distances, 95)),
        }

        print(f"Posterior predictive check: "
              f"median distance = {np.median(distances):.3f}, "
              f"p95 = {np.percentile(distances, 95):.3f}")

        return distances

    def counterfactual(self, intervention_fn, intervention_name="",
                       n_draws=50, n_replicates=5, seed=42):
        """Run a counterfactual experiment."""
        if self.posterior_samples is None:
            raise ValueError("Call calibrate() first")

        rng = np.random.default_rng(seed)
        indices = rng.choice(
            len(self.posterior_samples), min(n_draws, len(self.posterior_samples)),
            replace=False,
        )

        effects = []
        for idx in indices:
            theta = self.posterior_samples[idx]
            params_dict = dict(zip(self.param_names, theta))

            for _ in range(n_replicates):
                run_seed = rng.integers(0, 2**31)

                # Baseline
                base_data = self.simulator_fn(params_dict, seed=run_seed)
                base_summ = self.summary_fn(base_data)

                # Intervention
                mod_dict = intervention_fn(params_dict.copy())
                mod_data = self.simulator_fn(mod_dict, seed=run_seed)
                mod_summ = self.summary_fn(mod_data)

                effects.append(mod_summ - base_summ)

        effects = np.array(effects)
        self.provenance.setdefault("counterfactuals", []).append({
            "name": intervention_name,
            "n_draws": n_draws,
            "n_replicates": n_replicates,
            "mean_effect": effects.mean(axis=0).tolist(),
        })

        return effects

    def save_provenance(self, path):
        """Save full provenance record to JSON."""
        with open(path, "w") as f:
            json.dump(self.provenance, f, indent=2, default=str)
        print(f"Provenance saved to {path}")


# Assemble the pipeline
pipeline = CalibratedSimulatorPipeline(
    simulator_fn=lambda params, seed=42: toggle_switch_gillespie(
        ToggleSwitchParams(**params), t_max=200.0, seed=seed
    ),
    summary_fn=toggle_switch_summaries,
    param_names=["alpha_a", "alpha_b", "delta_a", "delta_b", "n"],
    prior_ranges={
        "alpha_a": (10, 100),
        "alpha_b": (10, 100),
        "delta_a": (0.1, 3.0),
        "delta_b": (0.1, 3.0),
        "n": (1.0, 5.0),
    },
)

# The full workflow:
# pipeline.set_observed_data(observed)
# particles, weights = pipeline.calibrate(n_populations=8)
# distances = pipeline.validate()
# effects = pipeline.counterfactual(double_delta_a, "double_delta_a")
# pipeline.save_provenance("toggle_switch_provenance.json")

The pipeline class encapsulates the science, but a standalone script is not yet a component that other tools can query; connecting it to the broader Discovery Workbench turns the calibrated simulator into an on-demand service for automated hypothesis testing.

9. Integrating with the Discovery Workbench

The calibrated simulator pipeline plugs into the Discovery Workbench as a component that accepts hypotheses (as parameter values or intervention specifications) and returns predictions with uncertainty estimates. This integration enables automated hypothesis testing loops where the Workbench generates candidate mechanisms, the simulator tests them, and the results guide the next round of hypothesis generation.

from discovery_workbench import Component, register


@register("calibrated_simulator")
class CalibratedSimulatorComponent(Component):
    """Discovery Workbench component for calibrated simulation.

    Exposes three actions:
    - predict: simulate from the posterior and return predictions
    - intervene: run a counterfactual experiment
    - compare_models: use ABC model selection to rank mechanisms
    """

    def __init__(self, pipeline):
        self.pipeline = pipeline

    def predict(self, n_samples=100, seed=42):
        """Generate posterior predictive samples."""
        rng = np.random.default_rng(seed)
        predictions = []

        for _ in range(n_samples):
            idx = rng.choice(len(self.pipeline.posterior_samples))
            theta = self.pipeline.posterior_samples[idx]
            params = dict(zip(self.pipeline.param_names, theta))
            result = self.pipeline.simulator_fn(
                params, seed=rng.integers(0, 2**31)
            )
            predictions.append(result)

        return predictions

    def intervene(self, intervention_fn, name=""):
        """Run a counterfactual and return effect with uncertainty."""
        effects = self.pipeline.counterfactual(
            intervention_fn, intervention_name=name
        )
        return {
            "mean_effect": effects.mean(axis=0).tolist(),
            "std_effect": effects.std(axis=0).tolist(),
            "ci_lower": np.percentile(effects, 2.5, axis=0).tolist(),
            "ci_upper": np.percentile(effects, 97.5, axis=0).tolist(),
            "intervention_name": name,
        }

    def provenance(self):
        return self.pipeline.provenance

10. Summary

This section demonstrated the complete arc of simulation-based scientific discovery. The pipeline built a stochastic simulator of a genetic toggle switch, generated synthetic experimental data, and designed summary statistics that capture bistable dynamics. ABC-SMC calibrated the simulator against the data. Posterior predictive checks validated the calibration, and counterfactual experiments tested mechanistic hypotheses. The pipeline class encapsulates this workflow with full provenance tracking, and the Discovery Workbench integration makes it available as a component in automated discovery loops. The same template applies to any stochastic simulator: swap the reaction network, the summary statistics, and the prior, and the calibration machinery transfers with minimal modification (though the choice of summary statistics and prior ranges requires domain-specific tuning for each new system).