Observation Models

Measure a simulation the way an experiment would: BOLD, EEG, functional connectivity and power spectra, declared as observation models.

Part of the running example, where stage 2 adds an observations block that declares what to measure.

Observation models define what to measure from a simulation: BOLD signals, functional connectivity, power spectral density and custom pipelines. Everything is defined in YAML: functions are declared in the functions: section and referenced by name in the observation pipeline:.

BOLD Signal

The Balloon-Windkessel hemodynamic model transforms neural activity into a BOLD fMRI signal. The pipeline is defined entirely in YAML, and each function is specified with its equation, callable, or time range, then chained together:

from tvbo import SimulationExperiment

exp_yaml = """\
dynamics:
    name: ReducedWongWang
    parameters:
        a: {value: 0.27}
        b: {value: 0.108}
        d: {value: 0.154}
        gamma: {value: 0.641}
        tau_s: {value: 100}
        w: {value: 0.6}
        J_N: {value: 0.2609}
        I_o: {value: 0.3}
    state_variables:
        S:
            equation: {rhs: "-S/tau_s + (1 - S) * gamma * d * (a*J_N*w*S + I_o + c_in) / (1 + exp(-1*(a*(J_N*w*S + I_o + c_in) - b)/d))"}
            initial_value: 0.5
    coupling_inputs:
        c_in: {}
integration:
    step_size: 0.1
    duration: 50000
functions:
    hrf_kernel:
        description: Hemodynamic response function (Volterra kernel)
        time_range:
            lo: 0
            hi: duration
            step: input.sample_period
        equation:
            rhs: "(1/3)*exp(-0.5*(t/1000)/tau_s)*sin(sqrt(1/tau_f - 1/(4*tau_s**2))*(t/1000))/sqrt(1/tau_f - 1/(4*tau_s**2))"
            parameters:
                tau_s: {value: 0.8}
                tau_f: {value: 0.4}
        arguments:
            - {name: duration, value: 20000.0}
        output: hrf_kernel_ts
    volterra_transform:
        description: Volterra nonlinear BOLD transformation
        equation:
            rhs: "(X - 1.0) * k_1 * V_0"
            parameters:
                k_1: {value: 5.6}
                V_0: {value: 0.02}
        output: bold_volterra
    subsample_bold:
        description: Subsample to BOLD TR
        equation:
            rhs: "subsample(X, stepsize)"
        arguments:
            - {name: stepsize, value: 10000}
        apply_on_dimension: time
        output: bold_subsampled
observations:
    - name: bold
      source: S
      period: 2000
      pipeline:
          - function: hrf_kernel
          - callable:
                module: scipy.signal
                name: fftconvolve
            output: bold_conv
            arguments:
                - {name: in2, value: hrf_kernel_ts}
                - {name: mode, value: full}
          - function: volterra_transform
          - function: subsample_bold
"""

exp = SimulationExperiment.from_string(exp_yaml)
print(exp.render_code('jax')[-1500:])
 on stderr, controlled by TVBO_LOG_LEVEL (default INFO).
    configure_logging()

    _parser = argparse.ArgumentParser(description="Run JAX-generated TVBO simulation")
    _parser.add_argument(
        "--spec",
        type=_Path,
        default=None,
        help="YAML experiment spec (default: ../spec/*.yaml next to this script)",
    )
    _parser.add_argument(
        "-o",
        "--output",
        type=_Path,
        default=None,
        help="Output directory for the result",
    )
    _args = _parser.parse_args()

    _spec = _args.spec
    if _spec is None:
        _candidates = sorted(
            (_Path(__file__).resolve().parent.parent / "spec").glob("*.yaml")
        )
        if not _candidates:
            raise SystemExit("No spec found; pass --spec PATH")
        _spec = _candidates[0]

    from tvbo import SimulationExperiment

    _experiment = SimulationExperiment.from_yaml(str(_spec))
    _state = _experiment.collect_state()
    _result = kernel(_state)
    logger.info(
        "Done: %s, shape=%s", type(_result).__name__, getattr(_result, "shape", None)
    )

    if _args.output is not None:
        _args.output.mkdir(parents=True, exist_ok=True)
        if hasattr(_result, "save"):
            _result.save(str(_args.output))
        else:
            import numpy as _np

            _np.savez(
                _args.output / "result.npz", data=getattr(_result, "data", _result)
            )
        logger.info("Wrote results to %s", _args.output)

The generated code shows each function definition followed by the pipeline composition. No Python functions are hardcoded: everything is generated from the YAML specification.

Note the emitted call subsample_bold(..., stepsize=20000). That count is not the value: 10000 written under arguments:, because sampling steps are resolved from the observation’s declared period (the BOLD TR, 2000 ms) divided by the integration step_size (0.1 ms). Declaring the TR in physical units keeps the pipeline correct when step_size changes; the literal under arguments: only supplies the function’s standalone default. Every backend routes through the same resolver, so the sample count is identical across jax, tvboptim, and tvb.

Pipeline Architecture

Observation pipelines chain functions sequentially. Each step’s output becomes the next step’s input automatically. Steps can be:

Type YAML Field Description
Equation equation.rhs Symbolic math expression, rendered to JAX
Callable callable: {module, name} Reference to an external library function
Kernel time_range + equation Generates a TimeSeries over a time range

Available Observation Models

Pre-defined observation models in tvbo/database/observation_models/:

Model File Description
BOLD (TVB) bold_tvb.yaml Balloon-Windkessel HRF with downsampling
FC FC.yaml Pearson correlation matrix

See Also