diff --git a/_posts/2025-10-25-phase-diagram-playbook.md b/_posts/2025-10-25-phase-diagram-playbook.md
new file mode 100644
index 0000000..d73c8d6
--- /dev/null
+++ b/_posts/2025-10-25-phase-diagram-playbook.md
@@ -0,0 +1,1675 @@
+---
+layout: distill
+title: "Shared Coordinates for Cross-Subject Brain Dynamics: Universal Latents and Directly Comparable Phase Diagrams"
+categories: [shared latent representations, tutorial, subject alignment, neural state-spaces, energy landscapes, intersubject comparability, state transitions, interpretable descriptors, brain dynamics]
+giscus_comments: true
+date: 2025-10-25
+featured: false
+
+authors:
+ - name: Julian Kędys
+ url: "https://www.linkedin.com/in/julian-kedys-a332222a6/"
+ affiliations:
+ name: Department of Digital Medicine, Poznań Supercomputing and Networking Center (PSNC), Polish Academy of Sciences
+ - name: Cezary Mazurek
+ url: "http://pl.linkedin.com/in/cezarymazurek"
+
+bibliography: 2025-10-25-phase-diagram-playbook.bib
+
+toc:
+ - name: "Overview"
+ - name: "Introduction"
+ subsections:
+ - name: "1. Big‑picture overview - why shared latents + PMEM → ELA → PDA?"
+ - name: "2. Where the methods come from (intuitive recap)"
+ - name: "3. Why shared latents are necessary (and useful)"
+ - name: "Motivation and contribution"
+ - name: "Pipeline overview"
+ subsections:
+ - name: "1) ELA-secure preprocessing (brief overview)"
+ - name: "2) Population-universal shared latent space"
+ - name: "3) Binarisation of latent time series"
+ - name: "4) Pairwise maximum-entropy (Ising) fitting"
+ - name: "5) Energy-Landscape Analysis (ELA): descriptors and kinetics"
+ - name: "6) Phase-Diagram Analysis (PDA): multi-observable placement"
+ - name: "Results in a word"
+ - name: "Robustness, uncertainty and diagnostics"
+ - name: "Reproducibility and artefacts"
+ - name: "Limitations and scope"
+ - name: "Outlook"
+
+---
+
+## Overview (TL;DR)
+
+We present a modular, modality-agnostic workflow that turns heterogeneous whole-brain time series into cohort-comparable, interpretable coordinates on a shared phase diagram, together with energy-landscape descriptors such as attractors, barriers, and kinetics.
+
+Key steps:
+1) population-universal latent spaces (Shared Response Model - SRM, Multiset Canonical Correlation Analysis - MCCA, Group PCA or Group ICA with consensus and automated dimensionality selection)
+2) per-latent binarisation to the +/-1 format
+3) Pairwise Maximum Entropy Model (PMEM) or Ising fitting: exact for small N, pseudo-likelihood, or variational Bayes
+4) energy landscape analysis (ELA): minima, disconnectivity, barriers, occupancies, kinetics
+5) phase diagram analysis (PDA): novel multi-observable placement on a shared reference surface with uncertainty
+
+Outputs include uncertainty, quality control, and interactive visuals. Methods are user-tweakable, reliable, and reproducible.
+
+---
+## Introduction
+
+#### 1. Big‑picture overview - why shared latents + PMEM → ELA → PDA?
+Modern whole‑brain recordings are heterogeneous across subjects, sessions, tasks and modalities. If we analyse each participant in their own idiosyncratic space, descriptors of “brain state” are not directly comparable. Our pipeline solves this in two moves:
+1. Population‑universal shared latents: We align subjects into a common, low‑dimensional space (SRM / MCCA / Group PCA‑ICA with consensus). Variables have stable meaning across participants and runs, so everything downstream is comparable and reproducible.
+2. Physics‑grounded descriptors: On the binarised latents we fit a pairwise maximum‑entropy model (PMEM/Ising), then read out two complementary summaries:
+
+ * Energy‑Landscape Analysis (ELA) - an attractor‑and‑barrier view of the fitted Ising energy. It yields minima/basins, disconnectivity graphs, barrier spectra, and kinetic descriptors (basin dwell times, mean first-passage times (MFPTs), committors, relaxation times). This is the mechanistic, state‑space view.
+ * Phase‑Diagram Analysis (PDA) - a macroscopic view that places each subject on the **($$\mu$$, $$\sigma$$)** plane of a disordered (via parametric perturbations) Ising model (SK‑like). In broad outline, it uses multiple observables at once to locate individuals relative to critical boundaries, providing cohort‑comparable coordinates and uncertainty.
+
+#### 2. Where the methods come from (intuitive recap)
+* **PMEM/Ising:** Among all binary models that match the empirical means and pairwise correlations, PMEM has maximum entropy. It is equivalent to the zero‑temperature Ising family used throughout statistical physics. Minimal assumptions; parameters are interpretable as fields $$h_i$$ and couplings $$J_{ij}$$.
+* **ELA:** Treat the fitted Ising as an energy landscape:
+
+ $$
+ E(\mathbf{s}) = -\sum_i h_i s_i - \tfrac{1}{2}\sum_{i\neq j} J_{ij} s_i s_j
+ $$,
+
+ over binary states $$\mathbf{s}\in\{-1,+1\}^N$$.
+
+ Local minima = attractors; energy differences = barriers; transition graphs + Markov kinetics = how the brain moves between them.
+* **PDA:** In spin‑glass models, the distribution of couplings matters. If the off‑diagonal $$J_{ij}$$ have mean $$\mu$$ and standard deviation $$\sigma$$ with $$h_i\approx 0$$, the system sits in regions (paramagnetic / spin‑glass / ferromagnetic) that govern ordering, glassiness and susceptibility. PDA maps each subject onto this phase surface so cohorts can be compared at a glance.
+
+#### 3. Why shared latents are necessary (and useful)
+* **Stable semantics:** Each latent represents the same population‑level pattern across participants, which makes ELA basins and PDA coordinates directly comparable.
+* **Tractability:** PMEM scales as $$\mathcal{O}(N^2)$$ parameters; a well‑chosen latent space puts us in the sweet spot between information retention and robust estimation.
+* **Downstream identifiability:** Binarisation → PMEM → ELA/PDA relies on on/off switching. Our alignment preserves this structure and gives us comparable switching rasters across the cohort.
+* **Utility‑forward:** with one aligned space we can publish shared phase diagrams and landscape reports that are re‑usable across datasets and modalities, enabling baseline‑independent comparisons and cross‑study synthesis.
+
+**Select practical advantages of our framework:**
+
+The workflow is largely standalone and implemented locally, from scratch - allowing researchers/analysts to adapt its workings to their exact needs - with numerical-stability and computational-efficiency improvements (relative to analogous implementations in the domain), novel algorithmic extensions (i.a. for multi-subject dimensionality reduction, comparative phase diagram analysis of heterogeneous brain dynamics), informative metrics and rich visualisations, and emphasis on **a)** automated parameter optimisation - not requiring domain-expertise or significant prior experience with the pipeline from the user, **b)** data-driven model selection, **c)** data-modality universality and independence of heuristics/meta-knowledge throughout the entire design process, and **d)** availability of alternative methods and hyperparameters for key processing/modelling stages, so as to best fit the needs of the user
+
+**Limitations worth remembering:**
+- Binarisation coarsens signals (whereas more granular discretisation becomes computationally prohibitive almost instantly for real-life data/problems)
+- Results depend on the selection of binarisation thresholds, dimensionality-reduction models and target counts of obtained latent features
+- For exact modelling methods, the set of possible states doubles in size with every additional feature/node/brain region
+- PDA assumes $$h\approx 0$$ and an SK‑like (Sherrington-Kirkpatrick) parametrisation
+- ELA explores single‑spin‑flip pathways, which is a limited and simplified assumption
+
+To counteract the influence of initial choice of (hyper)parameters, we quantify uncertainty in the workflow, track convergence wherever applicable, offer truly data-driven and bias-free optimisation of pipeline parameters, and expose diagnostics so these choices remain transparent and testable.
+
+---
+
+## Motivation and contribution
+
+Comparing whole-brain dynamics across individuals is hard without a common reference that preserves interpretability and quantifies uncertainty. This challenge becomes even more apparent in studies of complex brain processes spanning cognition, sensory integration, and perception; when the data are limited; or when comparing brain dynamics driven by systematically different sub-types of various neurodevelopmental, psychiatric, or neurodegenerative conditions. Aiming to address these challenges and reflecting the UniReps community's interest in representational alignment and comparability across subjects, datasets, and models, this post demonstrates:
+
+- a subject-alignment front-end that produces shared latent representations with stable semantics across a population, offering several alternative approaches
+- a stitched, physics-grounded back-end (PMEM to ELA to PDA) that yields mechanistically interpretable descriptors and shared phase-diagram coordinates derived with our original methodology
+- a robustness-first toolkit that includes custom-built consensus alignment, automated parameter selection, uncertainty quantification, diagnostics, and review-ready artefacts
+
+---
+
+## Pipeline overview
+
+1. Preprocess and engineer: ELA-secure detrending, conservative handling of missing data, safe within-subject run concatenation, per-region standardisation.
+2. Population-aware alignment and dimensionality reduction: shared latents via SRM, MCCA, or Group PCA or Group ICA with consensus; automatic dimensionality selection.
+3. Binarisation: per-latent threshold (usually median or mean, unless domain-expertise justifies, e.g., percentile-based thresholding), yielding +/-1 time series.
+4. PMEM or Ising fitting: exact (small N), safeguarded pseudo-likelihood, or variational Bayes - each enriched with its adequate set of solution-stability and significance/quality assessments.
+5. Energy-landscape analysis: attractors, barriers, disconnectivity graph, occupancies, kinetics, and many more descriptors providing mechanistic, biologically meaningful insight into brain dynamics, as well as facilitating direct and intuitive comparisons between subjects/cohorts.
+6. Phase-diagram analysis: multi-observable placement on a shared reference surface with our custom cost function, reports on confidence intervals, and near-criticality indices.
+
+**Data summary** (used for the development and testing of the pipeline; details not central to current discussion in themselves):
+
+Resting-state fUS from N=8 mice (7 Cre-lox ASD models spanning 4 subtypes; 1 control with no symptomatic manifestation modelled); 54 bilateral regions (27 L/R pairs; whole-brain collection) - unified across all the subjects; two runs per mouse (1494 frames for each recording session; TR ≈ 0.6 s); runs concatenated per subject.
+
+
+
+
+
+ Interactive 3D energy landscape for an example mouse
+
+
+
+
+
+
+
+
+ Interactive 3D energy landscape for another example mouse
+
+
+
+
+
+
+
+
+
+---
+
+## 1) ELA-secure preprocessing (brief overview)
+
+Aim: remove spikes, outliers, and slow global drifts while **preserving the on/off switching structure** that drives binarisation, PMEM fitting, and ELA/PDA. The procedure is modality-agnostic, non-invasive, and parameterised to be reproducible.
+
+* Adaptive per-region parameters are computed from simple statistics and Welch spectra, then re-adapted after each step if requested (robust mode).
+* Despiking uses derivative-based detection with local, percentile-scaled replacements in short contextual windows; consecutive spikes are handled as blocks.
+* Outlier removal is IQR-based with the same local, percentile-scaled replacement; an optional second pass is available.
+* Population-aware detrending uses a cohort-optimised LOESS trend (fraction selected by stationarity and autocorrelation reduction). The global trend is estimated on the mean signal and scaled per region, which corrects drift without flattening transitions.
+* Optional steps: light smoothing, bandpass filtering, breaking long flat runs, temporal standardisation, and spatial normalisation.
+* Outputs include per-step amplitude deltas and residual checks; we also report the concordance between pre- and post-detrending binary states to ensure switching patterns are retained.
+
+
+
+
+ Five-step preprocessing for a representative region (Striatum dorsal, R): Original → after despiking
+ (derivative-based detection with local replacement) → after outlier removal (IQR with local context) →
+ after universal LOESS detrending (global trend removed, transitions preserved) → after spatial
+ normalisation to [0, 1].
+
+
+
+
+
+
+
+ Cohort-wide evidence for a global trend. Top: absolute linear-trend correlation \(|r|\) for each
+ recording with a decision threshold (red dashed line). Bottom: relative trend magnitude with its
+ threshold. Many runs exceed both criteria, motivating a universal detrending step.
+
+
+
+
+
+
+ Model selection for detrending on an exemplar run. Top: global signal with linear, quadratic, and LOESS
+ fits (LOESS selected). Bottom: residuals after each method. LOESS yields the most stationary, least
+ autocorrelated residuals, hence chosen for the universal detrending step.
+
+
+
+
+
+
+ Region-wise example showing that global detrending preserves the on/off switching used downstream.
+ Left: raw vs detrended signals for dentate gyrus (top) and thalamus (bottom).
+ Right: binary rasters before/after median thresholding; “Concordance” in the panel titles reports the
+ fraction of timepoints whose binary state is unchanged by detrending.
+
+
+
+
+
+
+ Audit table across all recordings: linear-trend correlation and \(p\)-value, relative trend magnitude,
+ and stationarity of residuals after linear, quadratic, and LOESS detrending. Most datasets achieve
+ stationarity only with LOESS, supporting the universal-detrending choice.
+
+
+
+
+
+
+
+{% details Click to expand: practical notes %}
+
+* Replacement never injects artificial structure: values are drawn from local percentiles and clamped to local ranges; neighbors can be skipped to avoid bleeding flagged samples.
+* Block handling groups consecutive indices to avoid fragmentation; a de-blocking fix prevents long identical segments after replacements.
+* Universal LOESS fraction is chosen across the cohort to balance residual stationarity, autocorrelation reduction, and variance explained; region-wise application only scales that global trend.
+* All steps are seedable and config-driven; logs capture chosen parameters, max amplitude changes, and QC metrics for auditability.
+
+#### Remark:
+
+For ELA, methods must enhance quality without compromising temporal structure. Safe, beneficial steps:
+
+* Despiking: generally beneficial; percentile/MAD-based, no hardcoded amplitude thresholds.
+* Outlier removal: single or repeated; statistics-driven safeguards prevent over-removal.
+* Detrending: crucial for our data; prefer the universal LOESS approach for cross-subject/session/region comparability. Region-wise detrending should be used only with domain context; the cohort-wide analysis supports global detrending, with expert review welcomed.
+* Flat-blocks breaking: apply if replacements create long identical runs.
+* Smoothing: only very mild, primarily cosmetic; avoid aggressive windows that could distort binarisation.
+
+Bandpass filtering can be performed if not already applied to the data - many recording devices perform it automatically. Temporal standardisation and spatial normalisation are not required for ELA itself but are retained for general use; spatial normalisation is applied for the imputation pipeline to align signs and amplitudes. LOESS may shift some series below zero; this is expected and accounted for. All parameterisations are chosen to remain strictly ELA-compatible.
+
+ {% enddetails %}
+
+### Preprocessing Supplement I - biologically plausible imputation of time series for entire missing brain regions
+
+{% details Click to expand: Imputation %}
+
+
+### Brain-region imputation: why this sub-pipeline works
+
+This sub-pipeline fills in missing regional time series in whole-brain power-Doppler recordings while **preserving temporal structure** and providing **auditable quality checks**. It offers three complementary strategies and a common validation suite:
+
+### What the pipeline does (high-level summary)
+
+* **Detect & canonise**
+ Finds mice/runs with missing regions and re-indexes to a canonical (reference-mouse) ordering so every region has a stable label and posiion.
+
+* **Quantify bilateral coupling**
+ For each left/right pair it computes change-based correlation, raw correlation, directional agreement, magnitude coupling, lagged cross-correlation (with optimal lag), and frequency-domain coherence. These metrics tell us whether a contralateral trace is a **plausible shape donor** and provide thresholds for safe use.
+
+* **Offer three imputation routes**
+
+ 1. **Bilateral (shape-copy with lag alignment):**
+ Mirrors the contralateral region after shifting by a **median optimal lag** estimated from reference mice. It **does not scale amplitudes** (we work in normalised units), preserving the on/off **shape** that matters for downstream binarisation / PMEM / ELA/PDA. Optional light noise can be added in a non-deterministic mode.
+
+ 2. **Temporal (population median):**
+ Builds the missing series from the **median temporal pattern** of the same region across other mice (optionally across both runs). This is robust to outliers and yields stable reconstructions; with MAD-based jitter it can reflect natural variability while staying anchored to the cohort’s typical dynamics.
+
+ 3. **Clustering / nearest-neighbours:**
+ Groups reference mice/runs for the same region using correlation or Euclidean distance and imputes from the **cluster mean** of the nearest group. This conditions the reconstruction on **matched dynamics**, often outperforming global medians when cohorts are heterogeneous. A 3-D PCA visualisation makes the neighbourhoods and the imputed point **inspectable**.
+
+* **Validate, don’t guess**
+ Every imputed series is compared to the population using:
+
+ * mean absolute deviation to the cross-mouse mean,
+ * correlation with that mean,
+ * **coverage** within 1–2 s.d.,
+ * change-correlation, peak cross-correlation and optimal lag,
+ * mean coherence and peak-coherence frequency,
+ * **power-spectrum correlation**.
+ These metrics are printed and plotted, so each imputation carries its own **provenance and QC**.
+
+### Why this is useful for downstream physics-style analyses
+
+* **Shape-faithful**: methods preserve **temporal switching** structure (crucial for binarisation → PMEM → ELA/PDA).
+* **Cohort-aware**: temporal and clustering routes borrow information only from the **same labelled region** across other mice/runs.
+* **Bilateral advantage**: when hemispheric symmetry is strong, lag-aligned mirroring recovers realistic trajectories with minimal bias.
+* **Transparent & reproducible**: seeds are fixed; thresholds are explicit; NaNs and edge cases are handled defensively; outputs are re-indexed to a canonical order.
+* **Method flexibility**: you can pick the route that best matches your biology (e.g., bilateral for symmetric circuits; clustering for diverse cohorts) and still get the **same validation bundle**.
+
+---
+
+### Figures
+
+
+
+
+ Imputation outcome and bilateral synchrony (Postrhinal area, R; Run 1).
+ Top: the imputed trace (red) is plotted against reference mice (light blue) and the population mean (green); the y-axis is the normalised power-Doppler signal. Bottom: bilateral
+ change scatter (left minus right first differences) with the diagonal “perfect synchronisation” guide, a fitted trend line, and summary metrics (change-correlation, peak cross-correlation with its optimal lag, mean coherence). Together these panels show that the imputed series follows cohort-typical dynamics and remains temporally consistent with its contralateral partner.
+
+
+
+
+
+
+
+ Lag-structure and frequency-locking (Postrhinal area, R; Runs 1–2).
+ For each run, the left panel shows cross-correlation across lags (dashed zero line for reference), and the right panel shows magnitude-squared coherence versus frequency. Peaks indicate preferred temporal offsets and shared frequency content; stable peaks across runs support consistent reconstruction rather than overfit noise.
+
+
+
+
+
+
+
+ Validation summary.
+ For each run the pipeline reports deviation from the population mean, correlation with the mean, coverage within 1–2 standard deviations, change-correlation, peak cross-correlation and optimal lag, mean coherence, and power-spectrum correlation. These values provide an audit trail for every imputation.
+
+
+
+
+
+
+
+ Pipeline provenance.
+ The log records which regions are missing, cohort-level bilateral coupling benchmarks used as context, and the selected imputation route (here, clustering-based). This makes method choice and inputs explicit for later review.
+
+
+
+
+
+
+
+
+
+ Nearest-neighbours in PCA space (interactive).
+ Each marker is a reference mouse/run for the same region, embedded by PCA of the region’s time-series vectors; colours show clusters (nearest-neighbour groups). The red diamond is the imputed series; the blue marker(s) indicate the current mouse/run.
+
+
+
+{% enddetails %}
+
+
+### Preprocessing Supplement II — similarity assessment and conditional concatenation of recordings from different sessions for each mouse
+
+{% details Click to expand: Concatenation %}
+
+### Concatenating two recording runs per mouse
+
+This section documents the methodology we use to **decide whether two runs are similar enough to be concatenated**, how we **align** them when needed, and what diagnostics we plot to make the process auditable.
+
+### What “similar enough to concatenate” means
+
+Before concatenation, we compare the two runs **region-by-region** and turn the result into a single pass/fail decision.
+
+### 1) Regionwise similarity (time-domain patterns)
+
+For each region $$(r)$$ present in both runs, we check a scalar **similarity** between its run-1 series $$\mathbf{x}_r$$ and run-2 series $$\mathbf{y}_r$$. Several options are available:
+* **Pearson correlation**: requires $$\operatorname{corr}\big(\mathbf{x}_r,\mathbf{y}_r\big)\ \ge\ \tau_{\mathrm{corr}}.$$
+* **Spearman** (rank correlation): compute $$\rho\in[-1,1],$$ map to $$[0,1]$$ by
+ $$\rho_{01}=\tfrac{1}{2}(\rho+1),$$
+ and require $$\rho_{01}\ \ge\ \tau_{\mathrm{spearman}}.$$
+ *(Good for monotone but non-linear similarity)*
+* **Kendall’s $$\tau$$**: same mapping $$\tau_{01}=\tfrac{1}{2}(\tau+1),$$ require $$\tau_{01}\ \ge\ \tau_{\mathrm{kendall}}.$$
+* **Euclidean (shape) distance**: min–max normalise both series to $$[0,1]$$, compute
+ $$d=\frac{\lVert \mathbf{x}_r-\mathbf{y}_r\rVert_2}{\sqrt{T}},$$
+ and require $$d\ \le\ \tau_{\mathrm{dist}}.$$
+ *(Insensitive to absolute scale; tests waveform similarity)*
+* **Cosine similarity**: map $$\cos\theta\in[-1,1]$$ to $$[0,1]$$ via
+ $$c_{01}=\tfrac{1}{2}(\cos\theta+1),$$
+ and require $$c_{01}\ \ge\ \tau_{\mathrm{cos}}.$$
+
+Each region yields a boolean pass/fail; the resulting vector is our **regionwise mask**.
+
+### 2) (Optional) Distribution similarity
+
+As a distributional check, we run a **two-sample Kolmogorov–Smirnov test** per region and declare a pass when the $$p-value$$ exceeds a threshold, i.e. the two marginal distributions are not detectably different at the chosen level.
+
+### 3) Aggregating to a single gate
+
+We fuse regionwise pass/fail results into a final score using one of:
+* **fraction** (default): share of regions that pass; concatenate if $$\mathrm{fraction}\ \ge\ \tau_{\mathrm{agg}}.$$
+* **weighted**: same as above but weights each region by amplitude or variance.
+* **median / any / all**: robust/lenient/strict alternates.
+
+The gate is reported as:
+
+```
+[check_run_similarity] aggregator='', final_score=, pass=
+```
+
+Only if it **passes**, or if `force_concatenate=True`, do we proceed.
+
+### Alignment: making run-2 comparable to run-1
+
+If concatenation proceeds, we align **levels** of run-2 to run-1 **per region** (no temporal warping):
+
+* `alignment_mode='match_medians'` (default): for each region $$(r)$$,
+ $$
+ \mathbf{y}_r \leftarrow \mathbf{y}_r + \big(\operatorname{median}(\mathbf{x}_r)-\operatorname{median}(\mathbf{y}_r)\big)
+ $$
+* Alternatively, `match_means` uses the mean instead of the median.
+
+**Why this is safe:** the ELA/PDA pipeline is driven by **on/off switching and relative changes**. Shifting a time series by a constant to match central tendency **does not** distort the temporal structure we use downstream.
+
+### Preprocessing before the gate (optional but recommended)
+
+Both runs can first pass through a light, ELA-secure preprocessing stack (despiking, robust outlier handling, LOESS detrending, mild smoothing). Parameters are seedable and adapt across regions. This improves the reliability of similarity estimates without changing the switching dynamics that matter later.
+
+### Concatenation step
+
+After alignment, we intersect region indices and **horizontally concatenate** the two runs (time dimension doubles). An optional last **outlier smoothing** can be applied to the concatenated trace using a robust IQR rule.
+
+### Diagnostics and what the plots show
+
+The helper `show_intermediate_concat_plots(...)` produces a **3×2** panel (if preprocessing is enabled):
+
+* **Row 1:** Run-1 RAW (left) and Run-2 RAW (right). Orange dashed line = mean; bright-green dashed line = median.
+* **Row 2:** Run-1 Preprocessed (left) and Run-2 Preprocessed (right) with the same guides.
+* **Row 3:** **Aligned Run-2** (left; grey dashed = pre-alignment, green = aligned) and **Final Concatenated** (right; with mean/median guides).
+
+A shared legend clarifies the mean/median guides. These views make the gating, alignment, and final result fully inspectable.
+
+---
+
+### Figures
+
+
+
+
+ Run-similarity gate (console readout).
+ For each mouse, we compute a per-region similarity mask and aggregate it to a single decision.
+ The line shows the chosen aggregator, the final score, and pass/fail.
+ Only runs that pass are aligned and concatenated; failing pairs are skipped to avoid mixing incompatible sessions (although similarity checks can be manually overridden and the user could force concatenation anyway).
+
+
+
+
+
+
+ Concatenation diagnostics for a representative region.
+ Top row: raw run-1 (left, blue) and raw run-2 (right, red) with orange mean and bright-green median lines.
+ Middle row: the same region after preprocessing.
+ Bottom-left: alignment step—run-2 before (grey dashed) and after (green) median-matching to run-1.
+ Bottom-right: final concatenated trace with mean/median guides. These panels document that the two runs are comparable, that level alignment has worked as intended, and that the final time series is suitable for downstream analyses.
+
+
+
+---
+
+## Why this approach is robust and useful
+
+* **Prevents spurious mixing.** The gate stops concatenation when the two runs **do not** tell a consistent story for most regions. This protects subsequent ELA/PDA stages from artefactual discontinuities.
+* **Flexible similarity notions.** You can choose correlation-based (linear or rank), cosine (directional), or Euclidean-shape metrics, depending on whether absolute scale or monotonicity matters most for your data
+ The rank-based options (Spearman/Kendall) are especially **stable** for neural time series, where amplitude changes can be non-linear or non-stationary.
+* **Scale-safe alignment.** Median/mean level-matching fixes global offsets without altering temporal structure, keeping **on/off** switching intact—the key for binarisation and PMEM fitting.
+* **Transparent diagnostics.** The console summary and 3×2 figure make each decision inspectable. If a pair fails, you immediately see why; if it passes, you can verify that alignment has not introduced distortions.
+* **Configurable strictness.** Thresholds $$(\tau_{\mathrm{corr}},\ \tau_{\mathrm{spearman}},\ \tau_{\mathrm{dist}},\ \dots,\ \tau_{\mathrm{agg}})$$ and the aggregation rule control how strict the gate is. You can be conservative for clinical datasets and more permissive for exploratory work.
+
+---
+
+### Minimal recipe (what the functions do)
+
+1. **Preprocess** each run (optional, recommended).
+2. **Compute regionwise similarity** with your chosen method.
+3. **Aggregate** passes into a final gate score; stop if it fails.
+4. **Align levels** (`match_medians`/`match_means`) per region.
+5. **Concatenate** the two runs (time axis).
+6. **(Optional) Final outlier pass** on the concatenated series.
+7. **Plot diagnostics** for audit and reports.
+
+This procedure is **seeded, reproducible, and auditable**, and it keeps exactly the aspects of the signal that matter for your subsequent ELA/PDA pipeline.
+
+{% enddetails %}
+
+
+
+---
+## 2) Population-universal shared latent space
+
+Goal: obtain shared, ELA-secure latent components whose semantics are stable across subjects and runs, so downstream binarisation → PMEM → ELA/PDA is directly comparable and computationally tractable.
+
+#### Novelty factor at a glance:
+
+* Pick-and-mix reductions: (Group PCA, Group ICA, SRM, MCCA) with robust tweaks (varimax, multi-restart ICA, PCA-initialised SRM, whitened MCCA).
+* A population-aware, multi-metric objective (structural + temporal + method-specific) → auto-selection of dimensionality and hyperparameters.
+* Alignment toolkit (orthogonal Procrustes with column-wise sign fixes, Hungarian matching, and a neuro-Procrustes consensus).
+* Synergy-weighted consensus across methods with dominance checks, stability (RV) analysis, and per-component contribution audits .
+* Efficient, reproducible compute (SVD on averaged covariances, fast downsampling for temporal metrics, parallel grid search, seeded randomness).
+
+
+---
+### 2.1 Methods (population-aware, ELA-secure)
+
+All methods preserve temporal ordering (only linear projections or orthogonal transforms over channels), and we quantify temporal faithfulness (trustworthiness/continuity, autocorrelation preservation) .
+
+* Group PCA (variance capture + varimax): eigenvectors of the average subject covariance; optional varimax keeps components sparse/interpretable without breaking orthogonality .
+*Extras:* elbow detection, subject-wise EVR, reconstruction error, and post-rotation low-variance pruning.
+
+* Group ICA (shared independent subspaces): subject-wise PCA → concatenation → FastICA with multi-restart selection (best negentropy proxy via kurtosis deviation), independence verification (mean $$\lvert\mathrm{corr}\rvert$$, kNN mutual information , squared-corr checks), sign harmonisation of subject contributions before averaging.
+
+* SRM (shared timecourses with subject-specific maps): iterative orthogonal subject mappings $$(W_i)$$ and shared response $$(S)$$, PCA-based initialisation, Hungarian-matched alignment of mappings across subjects, orthogonality diagnostics, shared variance quantification .
+
+* MCCA (maximising cross-subject correlation): per-subject whitening → SVD on concatenated features (SUMCOR-style) → iterative refinement of subject loadings $$(a_i)$$ . We report cross-subject alignment, canonical correlations to shared response, orthogonality in whitened space, and shared variance in native space.
+
+
+
+
+
+ Group PCA scree. Light-blue bars show individual explained-variance ratio (EVR) per
+ component; the black line is cumulative EVR. The dotted vertical line marks the elbow (here at the 1st
+ component), and the title reports total variance explained (65.6%). This plot guides the initial range
+ for dimensionality selection before our multi-metric model choice.
+
+
+
+
+
+
+
+
+
+
+ SRM latent space.Left: consensus SRM mapping \(W\) (regions × components).
+ Colours indicate signed loadings (arbitrary overall sign but harmonised across subjects), highlighting
+ stable spatial patterns shared across the cohort. Right: example subject’s correlation matrix
+ between SRM latent time series. Light (near-white) off-diagonals indicate low cross-component correlation,
+ i.e. little superfluous overlap—evidence of efficient representation and dimensionality reduction.
+
+
+
+
+
+
+
+
+
+
+ MCCA latent space.Left: example subject’s correlation matrix between MCCA
+ components—again, light off-diagonals reflect low redundancy across latents. Right: subject-specific
+ MCCA projection matrix \(a\) (regions × components), showing how each brain region contributes to each
+ shared component; structured bands point to interpretable, population-aligned patterns.
+
+
+
+
+> Why four methods? They trade off interpretability, independence, correlation sharing, and variance capture. We score them with the same unified, ELA-aware metric suite and fuse them, yielding a stable, population-universal latent space.
+
+---
+
+### 2.2 Alignment & sign consistency (critical for comparability)
+
+* Orthogonal Procrustes with column-wise sign checks (flip a component if it anticorrelates with the reference) .
+* Hungarian matching aligns component order across methods/subjects by maximal absolute correlations .
+* Neuro-Procrustes consensus: iterative, order-robust alignment across methods (SVD-based reference) with final sign harmonisation.
+* Optional biological sign protocol (baseline-anchored flips) to stabilise polarities across datasets.
+
+---
+
+### 2.3 Metrics & selection (multi-objective, normalised to [0,1])
+
+Structural fidelity (RSA-style): correlation between original region-by-region correlation structure and that reconstructed from latents; sign preservation weighted by $$\lvert\mathrm{corr}\rvert$$; Procrustes disparity; Kruskal stress on inter-region distances.
+
+Temporal faithfulness: trustworthiness/continuity of neighbour relations, temporal neighbourhood hits, autocorrelation preservation at task-relevant lags.
+
+Method-specific:
+
+* PCA: mean subject EVR, reconstruction error.
+* ICA: independence (mean $$\lvert\mathrm{corr}\rvert$$↓, Mutual Information (MI)↓), kurtosis (≠3), sparsity.
+* SRM: shared alignment (Hungarian-matched), orthogonality, shared variance.
+* MCCA: cross-subject alignment, shared variance, canonical correlations, orthogonality (whitened).
+
+Normalisation uses principled maps (e.g., $$x \mapsto \tfrac{1}{1+x}$$ for “smaller-is-better”, linear $$[-1,1]\to[0,1]$$ for correlations).
+Composite score: weighted average (user- or default weights), then auto-optimise $$(k)$$ and hyperparameters via seeded, parallel grid search.
+
+{% details Click to expand: metric formulas %}
+
+Kruskal stress (upper-triangle, variance-matched):
+
+{% raw %}
+$$
+\mathrm{Stress}_1
+=\min_{b>0}\sqrt{\frac{\sum_{i
+
+
+ Cross-subject consistency of shared components. Each bar is the mean absolute correlation of the same
+ component across subjects (higher = more reproducible). Dashed/dotted lines show permutation and 95% CI thresholds.
+ Components above both guides generalise well across mice and are kept for the consensus latent space.
+
+
+
+
+
+
+ Alignment quality across methods. Bars show how well each consensus component matches its counterparts
+ from Group PCA, Group ICA, SRM and MCCA after Procrustes + Hungarian alignment (mean |corr|). Higher values and crossing
+ the dashed/dotted thresholds indicate robust cross-method agreement, supporting the fused consensus basis.
+
+
+
+
+
+
+
+ Method contributions to the consensus components. Stacked bars report the weighted influence of each
+ method (Group PCA / Group ICA / SRM / MCCA) on every consensus component. Balanced contributions mean the final space is
+ not dominated by a single method and retains structure that is consistent across approaches.
+
+
+
+
+
+
+
+ Example aligned loadings: Group PCA. Brain regions × components matrix (unit-normalised) after
+ cross-subject/method alignment. Warmer/colder cells mark stronger positive/negative loadings; light colours near zero
+ indicate sparse, non-overlapping contributions—useful for efficient representation and reduced redundancy.
+
+
+
+
+
+
+
+ Cross-method component correlations (post-alignment). Each block compares components across the four
+ methods. A sharp diagonal with light off-diagonals shows one-to-one matches and little superfluous overlap, validating
+ that different methods recover consistent latent directions.
+
+
+
+
+
+
+
+ Final consensus transformation. Regions × components loadings that define the population-universal
+ latent space used for all subjects. The mapping is sign-harmonised and unit-normalised so that latents have consistent
+ semantics across mice and runs.
+
+
+
+
+
+
+
+ Consensus latent space — within-subject orthogonality. Correlation matrix of consensus component
+ time courses for an example mouse. Near-zero off-diagonals (light colours) indicate low redundancy between latents,
+ which aids binarisation and stabilises the downstream PMEM/ELA/PDA steps.
+
+
+
+
+---
+
+{% details Click to expand: Mathematical underpinnings %}
+
+
+### 2.5 Mathematical cores (rigorous but compact)
+
+**SRM (orthogonal, correlation‑maximising form)**
+
+For subjects $$X_i \in \mathbb{R}^{T \times d}$$, find $$W_i \in \mathbb{R}^{d \times k}$$ with $$W_i^\top W_i = I$$ and a shared response $$S \in \mathbb{R}^{T \times k}$$:
+
+$$
+\min_{\{W_i\},\, S}\; \sum_i \lVert X_i W_i - S \rVert_F^2
+\;\Longleftrightarrow\;
+\max_{\{W_i\},\, S}\; \sum_i \operatorname{tr}\!\big(S^\top X_i W_i\big)
+\quad \text{s.t. } W_i^\top W_i = I .
+$$
+
+Updates:
+
+$$
+S \leftarrow \frac{1}{n} \sum_i X_i W_i ,
+$$
+
+then z‑score $$S$$ column‑wise.
+
+$$
+W_i \leftarrow U V^\top, \qquad \text{where } U \Sigma V^\top = \operatorname{SVD}\!\big(X_i^\top S\big).
+$$
+
+*The above formulation reflects our shared-T assumption - since all the subjects from the dataset used for developing the pipeline had the exact same number of frames in their time series.
+
+---
+
+**MCCA (SUMCOR‑style, whitened)**
+
+Let $$X_i$$ be centred and whitened to $$\tilde{X}_i$$. Find $$a_i \in \mathbb{R}^{d \times k}$$ maximising total cross‑correlation:
+
+$$
+\max_{\{a_i\}} \; \sum_{i
+
+
+
+ Binarisation sanity check: activity fraction over time.
+ The blue trace shows, at each time step, the proportion of latents in the +1 state after per-latent
+ median thresholding. The grey dashed line marks 0.5; orange dotted lines show the observed range
+ (min ≈ 0.20, max ≈ 0.80). The red dashed line is the fitted linear trend
+ (+0.052 percentage points per 100 steps; total change ≈ 1.57 pp). Near-stationary behaviour centred
+ around 0.5 indicates balanced on/off usage and supports downstream PMEM fitting; large drifts would
+ flag thresholding issues or residual global trends.
+
+
+
+
+---
+
+## 4) Pairwise maximum-entropy (Ising) fitting
+We model each time‑point’s binary latent vector $$\mathbf{s}\in\{-1,+1\}^N$$ with the pairwise maximum‑entropy (PMEM/Ising) distribution:
+$$
+P(\mathbf{s}) \propto
+\exp\Big(\sum_{i} h_i s_i + \sum_{i.
+
+### Inference routes (complementary, scale‑aware):
+* **Exact (small $$N$$)**: Enumerate all $$2^N$$ states to obtain the exact log‑likelihood and moments for gold‑standard checks
+* **Pseudo‑likelihood (PL)**: Optimise the sum of node‑wise logistic conditionals with L2 penalties and a safeguarded Armijo line‑search; enforce symmetry $$J_{ij}=J_{ji}, J_{ii}=0$$ and use a relative‑norm stopping rule (scale‑free, robust).
+
+ **Local field:**
+
+ $$
+ f_i^{(t)} = h_i + \sum_{j\neq i} J_{ij}s_j^{(t)} .
+ $$
+
+ **Node-wise conditional and log-conditional:**
+
+ $$
+ p\left(s_i^{(t)}\mid \mathbf s_{\setminus i}^{(t)}\right)
+ = \frac{\exp\big(s_i^{(t)} f_i^{(t)}\big)}{2\cosh f_i^{(t)}},
+ \qquad
+ \log p\left(s_i^{(t)}\mid \mathbf s_{\setminus i}^{(t)}\right)
+ = s_i^{(t)} f_i^{(t)}-\log\big(2\cosh f_i^{(t)}\big).
+ $$
+
+ **PL objective with L2 penalties (optimise over $$(h)$$ and the upper-triangle $$\{J_{ij}\}_{i
+VB-Ising diagnostics. Iteration summary and posterior uncertainties for fields/couplings. Large \(|z|=\lvert\mu\rvert/\sigma\) edges are strongly supported; high coefficient-of-variation flags potentially unstable couplings. Use these readouts to prioritise robust interactions when interpreting circuit-level dependencies.
+
+
+
+
+
+
+ Variational-Bayes PMEM: uncertainty and fit checks.
+ Panels summarise one VB run:
+ (i) posterior uncertainty vs coupling magnitude;
+ (ii, v) heatmaps of posterior s.d. for couplings (diagonal masked; log-scale variant);
+ (iii) (negative) ELBO trajectory decreasing across iterations (convergence);
+ (iv) histogram of field uncertainties;
+ (vi) ARD precision spectrum indicating data-driven shrinkage;
+ (vii) model-vs-empirical state probabilities (log–log; closer to the diagonal is better);
+ (viii) histogram of coupling coefficient-of-variation \( \mathrm{CV}=\sigma/|\mu| \);
+ (ix) two fit-quality indices (moment-matching accuracy). Together these quantify parameter credibility and goodness-of-fit before ELA/PDA.
+
+
+
+
+
+
+
+
+ Energy–probability diagnostic.
+ Empirical pattern probabilities \( P_{\mathrm{emp}}(\sigma) \) vs shifted energies \( E(\sigma)-E_{\min} \).
+ An approximately linear trend in log-probability vs energy (dashed fit; slope and \( R^2 \) shown) is consistent with Boltzmann structure.
+ Points are coloured by basin; the circled marker denotes the global minimum. Marginal histograms summarise energy and count distributions.
+ Deviations at the extremes flag rare states and help identify outliers before landscape and phase-diagram analyses.
+
+
+
+### Fit quality & sanity
+
+We report (i) moment matching (means and pairwise correlations), (ii) multi‑information explained and KL‑reduction vs. independence, and (iii) empirical‑vs‑Boltzmann pattern agreement.
+
+For Monte‑Carlo checks we use multi‑chain sampling with **R̂**/effective sample size (ESS) diagnostics and estimate observables (magnetisation $$m$$, Edwards–Anderson $$q$$, spin‑glass and uniform susceptibilities $$\chi_{\mathrm{SG}},\chi_{\mathrm{Uni}}$$, specific heat $$C$$).
+
+
+
+Pseudo-likelihood fit and sampling check. The I2/IN ratio and correlation \(r\) quantify global match to empirical statistics; the PL error markedly improves the independent baseline. BM–Metropolis runtimes confirm that the fitted model is tractable for sampling-based validation and energy-landscape analyses.
+
+
+
+Null-based quality control. Thresholds for first-moment error \(m_1\) and configuration-error \(C\) are estimated from nulls; all shown mice pass. This guards against over-fitting and ensures that downstream landscape metrics (barriers, relaxation) rest on statistically credible fits.
+
+
+
+Bootstrap accuracy with 95% CIs. For each mouse we report resampled accuracies and confidence intervals, quantifying the stability of the fitted model against sampling noise — a practical measure of reliability for comparative neuroscience.
+
+
+### Implementation highlights
+
+PL uses mean‑field initialisation, symmetric updates, Armijo backtracking and a relative gradient test; VB stores posterior precisions and ELBO traces for convergence auditing. (See robustness notes in Section Robustness, uncertainty and diagnostics.)
+
+---
+
+## 5) Energy-Landscape Analysis (ELA): descriptors and kinetics
+
+Once $$(h,J)$$ are fitted, the Ising energy
+
+$$
+E(\mathbf{s})=-\sum_i h_i s_i-\tfrac{1}{2}\sum_{i\neq j}J_{ij} s_i s_j
+$$
+
+induces a rugged energy landscape over $$\{-1,+1\}^N$$.
+
+We compute:
+* **Attractors and basins:** Local minima are states whose single‑spin flips all increase energy. Every state is assigned to a basin by steepest‑descent (or best‑improving) paths. We summarise basins by occupancies and a disconnectivity graph .
+* **Barriers and disconnectivity:** The barrier between basins $$(\alpha,\alpha')$$ is the minimum, over all paths connecting them, of the maximum energy along the path; the disconnectivity graph visualises these heights. Denoting the (symmetrised) minimal saddle by $$\overline{E}_{\alpha\alpha'}$$,
+
+$$
+\overline{E}_{\alpha\alpha'} = \min_{\gamma: \alpha\to\alpha'}\ \max_{\mathbf{s}\in\gamma} E(\mathbf{s}).
+$$
+
+* We also estimate a depth threshold (null‑model percentile) to guard against spurious minima.
+
+
+
+
+
+Distribution of pairwise barrier heights. Histogram of inferred barrier heights \(\overline{E}_{\alpha,\alpha'}\) between metastable basins in the fitted landscape; the dashed line marks the sample median. Mechanistically, more negative/lower barriers indicate easier inter-basin switches, whereas rarer higher barriers point to protected transitions. The spread quantifies heterogeneity of switching difficulty across the mouse’s brain-state repertoire.
+
+
+
+
+
+**Kinetics (Markov view):** Build a single-spin-flip Metropolis chain with proposal “flip one spin uniformly” and transition probability:
+
+$$
+p(\mathbf{s}\to\mathbf{s}^{(i)}) = \frac{1}{N}\min\{1,\ \exp\big(E(\mathbf{s})-E(\mathbf{s}^{(i)})\big)\},
+$$
+
+where $$\mathbf{s}^{(i)}$$ is $$\mathbf{s}$$ with spin $$i$$ flipped. This yields a $$2^N\times 2^N$$ transition matrix $$P$$ (or a restriction to the visited subgraph).
+
+**From $$P$$ we derive:**
+
+* Stationary distribution $$\pi$$, dwell-time distributions, and basin occupancy.
+* Mean first-passage times (MFPT): from a set $$A$$ to $$B$$ via the fundamental matrix $$Z=(I-Q)^{-1}$$ of the transient block .
+* Committors $$q_{AB}$$ solving $$(I-Q),q=b$$, where $$b$$ collects transitions into $$B$$ .
+* Relaxation spectrum (mixing time-scales): non-unit eigenvalues $$\lambda_i(P)$$ with $$\tau_i=-1/\log\lvert\lambda_i\rvert$$, and the Kemeny constant (mean mixing time) $$K=\sum_{i\ge 2}\frac{1}{1-\lambda_i}$$.
+
+
+
+
+Dwell-time distribution (empirical). The heavy right tail on a log-PDF axis is consistent with near-geometric escape from basins. Longer dwells reflect stabilised neural configurations (metastability), while short dwells reflect rapid exploration; the slope encodes an effective leaving rate.
+
+
+
+
+Slow relaxation spectrum. Relaxation times \(tau_i\) (from the fitted Markov operator) quantify how quickly perturbations along each mode decay. A dominant \(tau_1\) and a gap to subsequent modes signal slow inter-basin exchange and long memory in the dynamics — a hallmark of metastable neural regimes.
+
+
+
+
+
+
+ Empirical committor \(q_i\) for \(A=[150]\to B=[1309]\).
+ The committor is the forward-commitment probability that a state first hits \(B\) before \(A\).
+ Mass near \(q_i\approx 0\) indicates most visited states lie closer to \(A\); states with \(q_i\approx 0.5\)
+ are transition-like and highlight probable dynamical bottlenecks in the brain’s state graph.
+
+
+
+
+**Read-outs:** (i) attractor maps (patterns + labels), (ii) disconnectivity graphs, (iii) barrier distributions, (iv) transition/reachability matrices (one-step and multi-step), and (v) kinetic summaries (MFPT heatmaps, committor fields, relaxation spectra, Kemeny constants). These quantify stability, switching propensity, and heterogeneity of access between states.
+
+Crucially, these mechanistic and interpretable descriptors and metrics provide an additional high-level framework for comparing brain dynamics across different individuals, or even cohorts with systematically divergent patterns of neural activity - a discrete and more intuitive alternative to classic means for unifying/juxtaposing representations in computational systems.
+
+
+Energy-Landscape Analysis (ELA) — summary descriptors. The composite figure illustrates the standard read-outs used downstream of the fitted Ising model: (A)Local-minimum patterns (binary states for each attractor); (B)3-D energy surface with labelled minima (white dots) and most-probable transition paths (white arrows); (C)Direct transition counts between minima (Metropolis single-flip kernel); (D)Disconnectivity graph showing barrier heights that separate basins; (E)Basin visit frequencies (empirical occupancy); (F)Basin sizes (number of micro-states per basin in state-space); (G)Direct/indirect transition counts summarising multi-step reachability. Deeper basins and higher barriers indicate more stable, harder-to-leave states; denser transition lanes point to preferred switching routes.
+
+
+Basin graph (alternative A: mosaics). Each panel corresponds to one attractor (State 1–18). Circles denote individual binary microstates assigned to that basin; circle colour encodes Ising energy (cooler = lower). This compact view shows how densely each basin occupies nearby configurations and highlights within-basin heterogeneity (broader colour spread ⇒ greater internal energy variance).
+
+ Basin graph (alternative B: directed neighbourhoods). For selected attractors (States 5–8), nodes are microstates (colour = energy) and arrows indicate admissible single-spin-flip moves under the Metropolis kernel. The layered, fan-shaped structure reflects typical downhill funnels into each attractor; sparse cross-links indicate rarer exits via saddles.
+
+
+3-D energy landscape. A continuous rendering of the fitted Ising energy with numbered minima (basins) and a white transition skeleton connecting them through low-saddle routes. Valleys (cool colours) are deep, stable basins; ridges quantify barrier heights that regulate switching.
+
+
+2-D energy landscape. The same landscape as a contour map. This top-down view makes it easy to read relative heights along paths and to spot alternative routes between basins (branch points near saddles). Together with the 3-D view, it provides complementary intuition about basin depth (3-D) and path geometry (2-D).
+
+
+Mean first-passage times (MFPT). Heatmap of expected steps to reach each target state from each start state (both sorted by basin index). The bright diagonal reflects near-zero self-passage; block structure and asymmetric bands reveal easy vs hard cross-basin routes. White gridlines mark basin boundaries; the colour bar is in steps.
+
+
+Dwell-time distributions per basin. For each attractor basin, the violin shows the full distribution of residence times (frames) from the single-spin-flip dynamics; points display individual visits. Wider violins and higher medians indicate kinetically stable basins; narrow shapes near 1–2 frames indicate transient basins.
+
+
+Basin-visit raster over time. Colour-coded stripe plot of the visited basin label across the recording. Long same-colour blocks correspond to sustained dwell periods; frequent colour changes indicate rapid switching. This readout complements MFPT and dwell-time summaries by exposing the temporal ordering of visits.
+
+
+
+Model-fit quality and global kinetics. Reported are fit accuracies, stationary entropy \(H(\pi)\) (spread of state use), slow \(\tau_i\), spectral gap \(\lambda_1-\lambda_2\) (mixing speed), and the Kemeny constant (mean hitting time averaged over targets). Together they summarise how well the model matches the data and how swiftly the brain’s state dynamics explore the landscape.
+
+
+
+Boltzmann rank plot (basin-coloured). Empirical probabilities versus energy rank (log-scale) assess Boltzmann-like ordering: high-probability states sit among the lowest-energy configurations. Basin colours reveal which attractors dominate the high-probability tail; the dashed fit summarises the overall decay.
+
+
+
+---
+## 6) Phase-Diagram Analysis (PDA): multi-observable placement
+
+**Goal:**
+Place every subject on a *shared* Sherrington–Kirkpatrick‑like $$(\mu,\sigma)$$ phase surface using *multiple* observables at once, with uncertainty, so that cohorts become directly comparable without needing a fixed “healthy baseline”. PDA sits downstream of our shared‑latent → binarisation → Ising (PMEM) fit, and is designed to be robust, auditable, and reproducible from end to end.
+
+
+
+
+ Interactive 3D phase diagram with robust multi-subject positioning - illustrated for C as the example metric
+
+
+
+
+
+---
+
+### 6.1 What PDA estimates (observables and phase coordinates)
+
+For a fitted Ising model on binary latents $$\mathbf{s}\in\{-1,+1\}^N$$ with fields $$h_i$$ and couplings $$J_{ij}$$, PDA uses macroscopic observables:
+
+- **Magnetisation**
+ $$m=\frac{1}{N}\sum_i \langle s_i\rangle$$
+
+- **Edwards–Anderson order**
+ $$q=\frac{1}{N}\sum_i \langle s_i\rangle^2$$
+
+- **Spin‑glass susceptibility** (using the eigenvalues $$\{\lambda_k\}$$ of the spin covariance)
+ $$\chi_{\mathrm{SG}}=\frac{1}{N}\sum_{k=1}^{N}\lambda_k^2$$
+
+- **Uniform susceptibility**
+ $$\chi_{\mathrm{Uni}}=\frac{1}{N}\sum_{i\neq j}\mathrm{Cov}(s_i,s_j)$$
+
+- **Specific heat** (from energy variance)
+ $$C=\frac{\langle E^2\rangle-\langle E\rangle^2}{N},\qquad
+ E(\mathbf{s})=-\sum_i h_i s_i-\frac{1}{2}\sum_{i\neq j}J_{ij}s_is_j$$
+
+In practice, we compute these from Monte‑Carlo samples of the fitted model (with automatic convergence checks). For quick diagnostics, the first four are also obtainable directly from the binarised data.
+
+---
+
+### 6.2 The reference phase surface $$(\mu,\sigma)\mapsto\{m,q,\chi_{\mathrm{SG}},\chi_{\mathrm{Uni}},C\}$$
+
+We construct a *reference* coupling matrix $$J_{\mathrm{ref}}$$ (either *pooled* across the cohort or *control*‑only), then generate a dense grid over target $$(\mu,\sigma)$$ by *affinely transforming the off‑diagonal entries* of $$J_{\mathrm{ref}}$$. Let $$\mu_{\text{old}}$$ and $$\sigma_{\text{old}}$$ be the mean and std of off‑diagonal entries of $$J_{\mathrm{ref}}$$. For a target $$(\mu,\sigma)$$:
+
+$$
+J_{ij}^{(\mu,\sigma)}=
+\begin{cases}
+\big(J_{ij}-\mu_{\text{old}}\big)\dfrac{\sigma}{\sigma_{\text{old}}+\varepsilon}+\mu, & i\neq j,\\[6pt]
+0, & i=j~,
+\end{cases}
+$$
+
+with all **fields zeroed** $$h_i\equiv 0$$ (diagonal remains 0) to recover the canonical spin‑glass phase structure. For each grid point we Monte‑Carlo sample the observables above and cache five surfaces $$\{\mathrm{m},\mathrm{q},\chi_{\mathrm{SG}},\chi_{\mathrm{Uni}},C\}$$.
+
+**Working vs display grids:** We build a high‑resolution *working* grid used for optimisation/placement and optionally a wider *display* grid for visuals. We automatically pick $$(\mu,\sigma)$$ limits by running quick PL fits per subject to estimate native $$(\hat\mu,\hat\sigma)$$, then expand by a safety factor; we avoid overly coarse grids (e.g., if $$\Delta\mu$$ or $$\Delta\sigma>0.01$$).
+
+**MC convergence safeguards:** We run multiple chains with increasing sweep budgets until *all* observables reach $$\hat R<1.05$$ and an effective sample size threshold, or a maximum sweep cap is hit (in which case a warning is issued).
+
+---
+
+### 6.3 Multi‑observable placement (cost projection with balanced weights)
+
+Given a subject’s *observed* $$(m_o,q_o,\chi^{o}_{\mathrm{SG}},\chi^{o}_{\mathrm{Uni}})$$ from their binarised latents, we **project** onto the reference surface by minimising a *variance‑balanced* squared error:
+
+$$
+\underset{\mu,\sigma}{\arg\min}\;\sum_{k\in\{\mathrm{m},\mathrm{q},\chi_{\mathrm{SG}},\chi_{\mathrm{Uni}}\}}
+w_k\big(O_k^{\text{obs}}-\widehat{O}_k(\mu,\sigma)\big)^2,
+$$
+
+where $$\widehat{O}_k(\mu,\sigma)$$ is obtained by regular‑grid interpolation of the precomputed surfaces and weights $$w_k$$ default to the inverse *range* of each surface (less sensitive than $$1/\mathrm{var}$$). Optimisation uses L‑BFGS‑B on the working grid bounds; we also provide an **iso‑curve** fallback (intersect level sets of $$\chi_{\mathrm{SG}}$$ and $$\chi_{\mathrm{Uni}}$$ with a simple border‑safety check) and a brute‑force grid fallback.
+
+We return $$(\hat\mu,\hat\sigma)$$, the final cost value, and the method used (“cost_minimisation”, “iso_curves”, or “fallback_grid”).
+
+
+
+
+ Specific-heat landscape \(C(\sigma,\mu)\) with cohort placements.
+ The dashed ridge marks a near-critical band where responsiveness inflates. Subjects cluster on the gentler
+ slope below the ridge; the control mouse consistently occupies the lowest-\(\sigma\) (least heterogeneous
+ couplings), as expected. This panel anchors the multi-observable placement on a shared surface.
+
+
+
+
+
+
+
+ Cross-checks across observables.
+ A complementary 3-D view (top) with 2-D slices for q (left) and \(\chi_{\mathrm{SG}}\) (right).
+ Consistent subject ordering across panels indicates that placements are not driven by a single metric.
+ The control remains the lowest-\(\sigma\) point on the pooled reference.
+
+
+
+
+---
+
+### 6.4 Uncertainty, robustness, and diagnostics
+
+**Bootstrap CIs:** For each subject we run a *circular block bootstrap* along time, refit the Ising per resample, and recompute $$(\mu,\sigma)$$ or $$(m,q,\chi_{\mathrm{SG}},\chi_{\mathrm{Uni}})$$ as needed to report 95% intervals. Specific heat $$C$$ can be bootstrapped likewise via short MC runs per resample.
+
+**Control shrinkage (optional):** When a control is available, we stabilise its estimate by bootstrapping its $$J$$, pooling across the cohort, and forming a convex combination $$J^{\text{ctrl,shrunk}}=(1-\lambda)J^{\text{ctrl}}+\lambda J^{\text{pooled}}$$; we then summarise $$(\mu_c,\sigma_c)$$ from the shrunk $$J$$.
+
+**Cost bowls:** Around each optimum we display the *local cost landscape* (2‑D filled contour and 3‑D surface) to judge identifiability; narrow, well‑curved bowls suggest precise placement.
+
+**Group‑level tests:** Small helpers allow groupwise comparisons in $$\sigma$$ (e.g., Welch ANOVA and permutation tests) using the bootstrapped distributions.
+
+**Critical structure:** We plot *critical contours* (e.g., a fixed fraction of the maximum of an observable) on the display surfaces; a simple near‑criticality index is the minimal Euclidean distance from $$(\hat\mu,\hat\sigma)$$ to the chosen contour.
+
+
+
+
+
+ \(\mu\)–\(\sigma\) placements with uncertainty.
+ Dots show bootstrap means; ellipses are 95% confidence regions from circular block-bootstrap resamples.
+ Groupings separate along \(\sigma\) (heterogeneity), and the control consistently shows the
+ smallest \(\sigma\) on the pooled reference — a sanity check aligned with biological expectations.
+
+
+
+
+
+
+ Local 2-D cost surface (“identifiability map”).
+ Filled contours of the variance-balanced multi-observable discrepancy around the optimum in \((\sigma,\mu)\).
+ A tight, symmetric bowl indicates precise, well-conditioned placement; mild elongation reflects correlated
+ trade-offs between observables.
+
+
+
+
+
+
+
+ 3-D “cost bowl”.
+ The same objective shown as a 3-D surface. Clear curvature corroborates convergence and local identifiability.
+ (Together with the 2-D map above, this rules out flat minima or boundary locking.)
+
+
+
+
+---
+
+### 6.5 Practical recipe (what the code actually does)
+
+1) **Prepare data**: shared latents → binarise per latent to $$\pm1$$.
+2) **Rough bounds**: quick PL fits to get subject‑wise $$(\hat\mu,\hat\sigma)$$; expand to working ranges; verify resolution is fine enough.
+3) **Build reference surface**: choose `reference_mode ∈ {pooled, control}`, set $$h\equiv 0$$, sweep $$(\mu,\sigma)$$ by affine off‑diagonal transforms of $$J_{\mathrm{ref}}$$, run multi‑chain MC with convergence checks, and cache $$\{\mathrm{m},\mathrm{q},\chi_{\mathrm{SG}},\chi_{\mathrm{Uni}},C\}$$ surfaces.
+4) **Place subjects**: minimise the balanced multi‑observable cost (or use the iso‑curve / grid fallbacks with border guard).
+5) **Uncertainty**: bootstrap along time to report CIs; optionally shrink the control to summarise $$(\mu_c,\sigma_c)$$.
+6) **Visuals**: 2‑D contour panels and interactive 3‑D surfaces, with group‑colours and critical lines; “cost bowls” per subject for identifiability.
+
+---
+
+### 6.6 Mathematical and implementation notes (exact to this workflow)
+
+- **Affine mapping of $$J$$ to $$(\mu,\sigma)$$** modifies *only* off‑diagonals and preserves $$J_{ii}=0$$. This avoids artefactual self‑coupling and keeps the spin‑glass structure intact. Fields $$h$$ are explicitly zeroed for the surface.
+- **Observables from data vs model.** Fast “data‑side” $$m,q,\chi$$ provide immediate checks, while MC‑side estimates (including $$C$$) are used to build the surface; both routes are available.
+- **Balanced weights** $$w_k$$ default to inverse *range* (normalises sensitivity without over‑penalising heavy‑tailed metrics). Equal weights are also supported.
+- **Convergence gating** uses multi‑chain Gelman-Rubin $$\hat R$$ and a conservative ESS check with geometric back‑off of sweeps up to a cap; we report if the cap is reached.
+- **Border guard** prevents pathological “corner locking” when intersecting $$\chi$$‑iso‑curves on coarse grids.
+
+---
+### 6.7 Interpretation tips
+
+### Metric glossary:
+
+**$$\sigma$$ — network heterogeneity (dispersion of couplings)**
+**Meaning:** standard deviation of off‑diagonal $$J_{ij}$$ on the reference surface.
+**High:** uneven, subnetwork‑biased interactions; rugged energy landscape with many competing basins; dynamics readily reconfigured by *local* nudges.
+**Low:** broadly uniform interactions; smoother coordination; fewer competing minima (less glassy).
+
+**$$\mu$$ — net coupling / co‑activation balance (mean coupling)**
+**Meaning:** mean of off‑diagonal $$J_{ij}$$ (with $$h\approx 0$$ on the surface).
+**High (positive):** stronger global ordering/co‑activation; coherent whole‑system shifts are easy.
+**Low/negative:** interactions cancel or oppose; regions act more independently/segregate.
+
+**$$m$$ — magnetisation (whole‑brain activation bias)**
+**Definition:** $$m=\tfrac{1}{N}\sum_i \langle s_i\rangle$$.
+**Large $$|m|$$:** prolonged hypo‑ or hyper‑activation (tonic bias; long runs in one sign).
+**Near 0:** balanced on/off usage.
+**Note:** on the reference surface $$h=0$$, so any non‑zero $$m$$ reflects ordering from $$\mu$$ (or thresholding bias when computed directly from data).
+
+**$$q$$ — Edwards–Anderson order (pattern rigidity)**
+**Definition:** $$q=\tfrac{1}{N}\sum_i \langle s_i\rangle^2$$ (persistence per unit regardless of sign).
+**High:** recurring, “frozen” configurations; can be rigid even when $$m\approx 0$$ (symmetry‑related states).
+**Low:** flexible/exploratory dynamics with weak per‑unit bias.
+
+**$$\chi_{\mathrm{SG}}$$ — spin‑glass susceptibility (sensitivity to *local* perturbations)**
+**Meaning:** response to heterogeneous, small nudges; equals the sum of squared eigenvalues of the spin covariance (normalised by $$N$$).
+**High:** small local changes can reconfigure the network; hallmark of glassy, high‑$$\sigma$$ regimes with many shallow basins.
+**Low:** locally robust; topology resists piecemeal perturbations.
+
+**$$\chi_{\mathrm{Uni}}$$ — uniform susceptibility (sensitivity to a *global* nudge)**
+**Meaning:** response when all units are pushed equally (normalised sum of pairwise covariances).
+**High:** easy, coherent whole‑brain shift (often increases with $$\mu$$, decreases with $$\sigma$$).
+**Low:** globally inertial/decoupled.
+
+**$$C$$ — specific heat (breadth of accessible repertoire)**
+**Definition:** $$C=\big(\langle E^2\rangle-\langle E\rangle^2\big)/N$$, i.e., energy variance per unit.
+**High:** wide repertoire; peaks near phase boundaries/critical bands (heightened responsiveness).
+**Low:** narrow variability; stereotyped dynamics.
+
+---
+
+### Typical regimes (reading combinations):
+
+* **Ferromagnetic‑like:** $$\mu$$ high, $$\sigma$$ low → large $$\|m\|$$, high $$q$$, high $$\chi_{\mathrm{Uni}}$$, low $$\chi_{\mathrm{SG}}$$; $$C$$ can peak near the ordering boundary.
+* **Spin‑glass‑like:** $$\mu\approx 0$$, $$\sigma$$ high → $$m\approx 0$$ but $$q$$ elevated, $$\chi_{\mathrm{SG}}$$ high, $$\chi_{\mathrm{Uni}}$$ low; many metastable basins and irregular switching.
+* **Paramagnetic‑like:** $$\mu\approx 0$$, $$\sigma$$ low → $$m\approx 0$$, low $$q$$, both susceptibilities low, $$C$$ low; weak coordination, noise‑like exploration.
+* **Near‑critical band:** $$C$$ high with elevated susceptibilities; large fluctuations and long correlation lengths.
+
+---
+
+### Practical notes:
+
+* PDA surfaces set $$h\approx 0$$; placements therefore reflect coupling structure $$\big(\mu,\sigma\big)$$ rather than tonic biases.
+* Compare subjects on the **same** reference (pooled vs control can shift absolute positions). Use uncertainty ellipses/cost bowls to judge identifiability.
+* Large $$\|m\|$$ in data space can arise from binarisation thresholds or residual trends—validate preprocessing.
+* PDA gives macroscopic coordinates; pair with ELA for basin‑level mechanisms (attractors, barriers, kinetics).
+
+---
+
+
+
+
+ Multi-observable phase surfaces with subject placements.
+ Each panel shows an observable over \( (\sigma,\mu) \): top row — magnetisation \( m \), Edwards–Anderson order \( q \), spin-glass susceptibility \( \chi_{\mathrm{SG}} \); bottom row — uniform susceptibility \( \chi_{\mathrm{Uni}} \) and specific heat \( C \).
+ Points mark each subject’s placement on the pooled reference. The steep wing indicates a near-critical band; the control occupies the smallest \( \sigma \).
+
+
+
+
+---
+
+### 6.8 Reproducibility knobs (defaults)
+
+- MC: chains = 12, start sweeps = 8k, burn‑in = 1k, cap = 128k, $$\hat R$$ tol = 1.05.
+- Working grid: $$140\times140$$ by default; guard: $$\max(\Delta\mu,\Delta\sigma)\le 0.01$$.
+- Objective: weights = “balanced” (inverse range), optimiser = L‑BFGS‑B within bounds.
+- Bootstrap: circular blocks, user‑set size; control shrinkage parameter $$\lambda\in[0,1]$$.
+
+---
+
+### 6.9 Illustrative pseudocode (orientation only — code is not released here)
+
+> The snippet below mirrors the sequence used to generate the figures; it is descriptive rather than an API.
+
+```python
+# 1) Build the reference phase surfaces (pooled or control)
+phase = build_phase_surfaces(
+ binaries_pm1 = BINARIES_PM1_OR_01, # dict: subject → DataFrame (±1)
+ reference = {"mode": "pooled", "subject": None}, # or {"mode": "control","subject": "ID"}
+ grid_size = 150, # working grid resolution
+ mc_settings = {"chains": 12, "burn_in": 1000, "start_sweeps": 8000}
+)
+
+# 2) Place subjects via multi-observable cost (balanced weights by default)
+positions = place_subjects_on_surface(
+ binaries_pm1 = BINARIES_PM1_OR_01,
+ phase = phase,
+ method = "cost", # or "iso_curves"
+ weights_mode = "balanced"
+)
+
+# 3) Inspect local identifiability (“cost bowl”) for one subject
+w = objective_weights(phase, mode="balanced")
+plot_cost_landscape(
+ subject_id = "SUBJECT_ID",
+ positions = positions,
+ phase = phase,
+ window = 0.06,
+ steps = 200,
+ weights = w
+)
+
+```
+
+---
+
+### 6.10 Caveats specific to PDA
+
+- PDA assumes an SK‑like parameterisation (coupling distribution matters); it analyses **macrostates** and does not replace ELA’s mechanistic, basin‑level descriptors.
+- The choice of reference surface (pooled vs control) can shift placements; we therefore expose the cost bowls and allow both options to be reported.
+- Grid resolution and MC budgets matter near sharp boundaries; guards and diagnostics make this explicit.
+
+---
+
+## Results at a glance
+On resting-state functional ultrasound (fUS) recordings (mesoscopic, whole-brain), we observe low placement residuals $$(10^{-6}–10^{-4})$$, tight bootstrap confidence regions, convergent models, and stable ordering under pooled vs subgroup phase references; example estimates span **σ ≈ 0.15–0.32** and **μ ≈ −0.01 to +0.03**. The outputs - susceptibility to perturbations, ordering vs glassiness, transition propensity - form compact, biologically meaningful fingerprints.
+
+For experimentalists, this is a mechanistic dashboard (what
+states exist, how deep, how likely to switch); for theorists, it anchors subjects in a physics-grounded phase
+space with interpretable axes.
+
+Overall, the multi-observable placement and the combination of shared embeddings with ELA/PDA provide reliability-minded, comparable, and interpretable read-outs that support discovery, phenotyping, and model-based hypothesis generation across cohorts, data modalities, tasks, and species.
+
+---
+
+## Robustness, uncertainty and diagnostics
+
+- Uncertainty via variational-Bayes posteriors for h and J; bootstrap intervals for mu and sigma and for near-criticality; block bootstrap for autocorrelation
+- Convergence checks including MCMC R̂ and effective sample size where applicable, pseudo-likelihood relative-norm stopping, and ELBO improvements for variational Bayes
+- Quality-control gates including permutation or null thresholds for spurious minima, component-stability filters, and sensitivity to number of latents and alignment choice
+- Ablations: pooled versus subgroup reference surfaces; method toggles among SRM, MCCA, and ICA; median versus mean thresholds
+
+---
+## Reproducibility and artefacts
+
+We report key settings and diagnostics to make the computational provenance clear, even though the code is not released with this post. Runs are seed-controlled with machine-parsable configuration files; convergence (e.g., R̂, ESS) and grid resolution checks are documented in the figure captions and text.
+
+**Example configuration stub (indicative):**
+
+```yaml
+# example-config.yaml
+seed: 123
+preprocess:
+ detrend: conservative
+ concat_runs: true
+ impute: short_gaps_only
+alignment:
+ methods: [SRM, MCCA, GroupPCA, GroupICA]
+ consensus: true
+ select_dim: auto
+binarise:
+ threshold: median
+ising:
+ mode: PL # {EXACT|PL|VB}
+ l2_h: 1e-5
+ l2_J: 1e-4
+ pl_tol: 1e-6 # relative-norm stopping
+ vb:
+ prior_prec_h: 6.0
+ prior_prec_J: 30.0
+ela:
+ minima_search: exhaustive
+ kinetics: true
+pda:
+ observables: [m, q, chiSG, chiUni, C]
+ bootstrap: true
+ ci_level: 0.95
+reports:
+ phase_report: true
+ landscape_report: true
+ ```
+
+---
+
+## Limitations and scope
+
+- Binarisation coarsens signals, but enables interpretable PMEM fitting and stable cross-subject comparability
+- Latent selection is influential; consensus and metric-guided selection mitigate but semantics remain partly model-dependent
+- The choice of reference surface affects PDA; we quantify sensitivity and expose cost bowls for transparency
+- Designed for resting or task neural time series; extension to other binarisable dynamical systems is often straightforward
+
+---
+
+## Outlook
+
+Population-universal latents combined with physics-grounded descriptors provide a shared language for multi-subject brain dynamics that is portable across modalities, tasks, and species, and a bridge to mechanistic interpretation and clinical translation. Planned extensions include multi-modal fusion, alignment-aware causal probes, truly dynamic/directional extensions of the methodology (e.g., by incorporating Langevin-based methods and attractor networks), developing modular workflows for modelling the state-spaces of consecutive processing stages in the brain under cognitive tasks, and targeted clinical studies.
+
+---
+
+## Appendix: Mathematical details
+
+Below are the core mathematical underpinnings of the pseudo‑likelihood and variational‑Bayes alternatives to the exact‑likelihood pairwise maximum‑entropy model (PMEM / Ising), as well as of the relevant auxiliary methods (e.g., for VB, the ARD/shrinkage or convergence diagnostics).
+
+{% details Click to expand the entire Appendix %}
+
+## Notation (common to PL and VB)
+
+- Binary state $$\mathbf{s}\in\{-1,+1\}^N$$ with samples $$\{\mathbf{s}^{(t)}\}_{t=1}^T$$.
+- Parameters $$\theta = \big[h_1,\ldots,h_N,\ \{J_{ij}\}_{i
+
\ No newline at end of file
diff --git a/assets/plotly/2025-10-25-phase-diagram-playbook/3D-ELA2.html b/assets/plotly/2025-10-25-phase-diagram-playbook/3D-ELA2.html
new file mode 100644
index 0000000..5e836b5
--- /dev/null
+++ b/assets/plotly/2025-10-25-phase-diagram-playbook/3D-ELA2.html
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/assets/plotly/2025-10-25-phase-diagram-playbook/3D-PDA.html b/assets/plotly/2025-10-25-phase-diagram-playbook/3D-PDA.html
new file mode 100644
index 0000000..e087cd3
--- /dev/null
+++ b/assets/plotly/2025-10-25-phase-diagram-playbook/3D-PDA.html
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/assets/plotly/2025-10-25-phase-diagram-playbook/3D-brain-regions-Louvain.html b/assets/plotly/2025-10-25-phase-diagram-playbook/3D-brain-regions-Louvain.html
new file mode 100644
index 0000000..5bcdf4e
--- /dev/null
+++ b/assets/plotly/2025-10-25-phase-diagram-playbook/3D-brain-regions-Louvain.html
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/assets/plotly/2025-10-25-phase-diagram-playbook/3D-impu.html b/assets/plotly/2025-10-25-phase-diagram-playbook/3D-impu.html
new file mode 100644
index 0000000..b52759b
--- /dev/null
+++ b/assets/plotly/2025-10-25-phase-diagram-playbook/3D-impu.html
@@ -0,0 +1,2 @@
+