TDE-RICA: Time-Delay Embedded Reconstruction Independent Component Analysis

Extract shared spatiotemporal motifs from multi-sample whole-brain neural activity recordings.

Based on Toyoshima et al. (2024), "Ensemble dynamics and information flow deduction from whole-brain imaging data", PLOS Computational Biology 20(3): e1011848. doi:10.1371/journal.pcbi.1011848

Overview

TDE-RICA combines Time-Delay Embedding (TDE) with Reconstruction Independent Component Analysis (RICA) to decompose neural activity time series into a small set of interpretable spatiotemporal motifs and their time-varying occurrence weights.

Key Capabilities:

Extract shared dynamical motifs across multiple samples (animals)
Handle missing data via matrix factorization imputation
Cross-validate model hyperparameters (embedding dimension, number of components)
Compare real vs simulated neural dynamics
Comprehensive time-series similarity analysis (alignment, distribution, continuous dynamics)

Input: [T x N x S] tensor of normalized neural activity
Output: Motifs [num_comp x dim_embed x N] and occurrences [T_embed x S x num_comp]

Repository Structure

TDE-RICA/
├── README.md                          # This document
├── LICENSE                            # MIT License
├── setup.py                           # Python package setup
│
├── +tderica/                          # MATLAB package (refactored)
│   ├── delayembed.m                   # Time-delay embedding
│   ├── delayembed_inv.m              # Inverse TDE
│   ├── select_subset.m               # Clean subset selection
│   ├── decompose.m                   # TDE-RICA core
│   ├── matrix_factorize.m            # L-BFGS imputation
│   ├── sort_cluster.m                # Hierarchical clustering
│   ├── visualize.m                   # Figure generation
│   ├── project.m                     # Component projection
│   └── compare_sim_to_real.m         # Real vs sim comparison
│
├── tderica/                           # Python package
│   ├── __init__.py                    # TdeRICA class + exports
│   ├── _tde.py                        # TDE + FNN dimension estimation
│   ├── _select_subset.py              # Greedy NaN-free selection
│   ├── _rica.py                       # RICA + FastICA fallback
│   ├── _decompose.py                  # Core decomposition pipeline
│   ├── _matrix_factorize.py           # L-BFGS + ALS solvers
│   ├── _sort_cluster.py               # Clustering + optimal leaf ordering
│   ├── _project.py                    # Least-squares projection
│   ├── _compare.py                    # Real vs sim scatter matrix
│   ├── _validate.py                   # Cross-validation tools
│   ├── _analysis.py                   # Occurrence similarity analysis
│   ├── _similarity.py                 # Comprehensive metrics suite
│   └── _visualize.py                  # Python visualization
│
├── tderica_main.m                     # MATLAB entry point (current)
├── tderica_test.py                    # Python unit tests (14 tests)
├── test_matlab_alignment.py           # MATLAB-Python alignment tests (32 tests)
├── tderica_orders.mat                 # Pre-computed cluster orders
│
└── demo_compare_real_and_simulation.m # Demo script

Installation

Python Package

Requires Python >= 3.10, NumPy >= 1.23, SciPy >= 1.7.

pip install -e .

Optional dependencies:

matplotlib -- visualization
scikit-learn -- FastICA fallback
joblib -- parallel TDE

MATLAB Toolbox

Add the repository root to your MATLAB path:

addpath('/path/to/TDE-RICA');

Requires MATLAB R2019b+ with Statistics and Machine Learning Toolbox (rica, linkage, optimalleaforder).

GNU Octave Compatibility

The MATLAB code also runs under GNU Octave (tested on 11.3.0) via a small compatibility layer. Octave lacks several MATLAB-only features used by TDE-RICA — rica, the internal classreg.learning.fsutils.Solver, histcounts/histcounts2, xcorr, alpha, and MATLAB-shaped plotmatrix outputs. octave_compat/ provides Octave replacements and setup_octave.m wires everything up:

% once per Octave session (from the repository root)
setup_octave

% then run the MATLAB entry point exactly as in MATLAB
tderica_main(tcrsNArranged, strNames, 300, 14, 1);

Prerequisites — the Octave statistics package (provides pdist, squareform, linkage, optimalleaforder):

octave --eval "pkg install -forge datatypes"
octave --eval "pkg install -forge statistics"

What the compatibility layer provides

MATLAB-only feature	Octave equivalent
`rica(...)`, `transform(mdl,X)`, `mdl.TransformWeights`	`octave_compat/rica.m` — Reconstruction ICA (Le et al. 2011 objective, L-BFGS with analytic gradient). The statistics package also ships a `rica`; `setup_octave.m` puts `octave_compat/` first so TDE-RICA's exact call (`'NonGaussianityIndicator'`) is accepted.
`classreg.learning.fsutils.Solver` (L-BFGS)	`octave_compat/+classreg/+learning/+fsutils/Solver.m` backed by `octave_compat/lbfgs.m`
`histcounts`, `histcounts2`	`octave_compat/histcounts*.m` (+ `binidx.m`)
`xcorr(...,'normalized',maxlag)`	`octave_compat/xcorr.m` — FFT-based, MATLAB-compatible lags and normalisation
`alpha(A)` on images	`octave_compat/alpha.m` — sets image `AlphaData`
MATLAB-shaped `plotmatrix` outputs (`AX(i,j)`, `Hax(i)`)	`octave_compat/plotmatrix.m` — shadows Octave's flat-vector variant
`get(H,'BinEdges')` from plotmatrix histograms	handled in `visualize.m` with a `histcounts` fallback (Octave has no `histogram` objects)

Verified differences / notes

The Octave RICA minimises 0.5*‖X − X·W·W'‖² + λ·Σ log cosh(X·W) (no whitening) so the reconstruction identity X ≈ occurrences·motifs holds exactly. On real data its reconstruction error equals the rank-K PCA baseline to 4 decimal places (e.g. 0.9262 vs 0.9262 for dim_embed=60, n_components=5).
The Python reference port's RICA (whiten=True default) was found to under-reconstruct on the same real data (error 1.10 > PCA 0.93); use whiten=False there if you need the Python baseline to match reconstruction quality.
.fig saving under Octave emits a harmless "unable to save onCleanup variables" warning; the file is still written.

Quick Start

Python API

import numpy as np
from tderica import TdeRICA

# data: [T x N x S] tensor (time x neurons x samples), may contain NaN
model = TdeRICA(dim_embed=300, n_components=14, lambda_=0.001)
model.fit(data)

occurrences = model.occurrences_   # (T_embed, S, num_comp)
motifs      = model.motifs_        # (num_comp, dim_embed, N)
reconstruction = model.reconstruct()  # (T, N, S)

MATLAB API

% tcrsNArranged: [T x N x S] normalized activity
% strNames:      {N x 1} cell array of neuron identifiers
tderica_main(tcrsNArranged, strNames, 300, 14, 1);

Algorithm Details

Core Pipeline

The TDE-RICA pipeline consists of four main stages:

graph TD
    A[Input: T x N x S tensor] --> B[Step 1: Select Clean Subset]
    B --> C[Step 2: Time-Delay Embedding]
    C --> D[Step 3: RICA Decomposition]
    D --> E[Step 4: Matrix Factorization Imputation]
    E --> F[Step 5: Sorting & Clustering]
    F --> G[Output: Motifs + Occurrences]

Step 1: Time-Delay Embedding (TDE)

Converts each neuron's instantaneous activity into a dim_embed-length trajectory window.

MATLAB:

for p = 1:dimEmbed
    tcrsEmbed(:, p, :) = tcrs(p:end-dimEmbed+p, :);
end

Python:

for p in range(dim_embed):
    out[:, p, :] = tcrs[p:T - dim_embed + 1 + p, :]

Inverse TDE reconstructs the original time series by averaging overlapping contributions (MATLAB: mean(..., 'omitnan'); Python: np.nanmean).

Step 2: Reconstruction ICA (RICA)

Learns independent spatiotemporal motifs by minimizing:

$$\min_{W,M} |X - WM|_F^2 + \lambda \sum_i \mathbb{E}[g(w_i^T x)]$$

where $g(\cdot)$ is a super-Gaussian contrast function (log-cosh), $W$ are occurrence weights, and $M$ are motif patterns.

Key difference from standard ICA: TDE-RICA assumes the weights of temporal patterns are independent (transposed setup), not the temporal patterns themselves.

MATLAB: Uses rica() from Statistics and Machine Learning Toolbox with NonGaussianityIndicator=ones(1,numComp) for super-Gaussian sources.

Python: Custom L-BFGS implementation with optional whitening, matching MATLAB's objective.

Step 3: Matrix Factorization for Imputation

Extends the clean-subset decomposition to all samples and cells by solving a block matrix factorization:

$$\min_{W_{add}, M_{add}} \left| \begin{bmatrix} W_{base} \ W_{add} \end{bmatrix} \begin{bmatrix} M_{base} & M_{add} \end{bmatrix} - \begin{bmatrix} X & A \ B & C \end{bmatrix} \right|_F^2$$

with $W_{base}$ and $M_{base}$ fixed.

Solvers:

L-BFGS (default): Joint optimization using curvature information
ALS (Python only): Alternating least squares with closed-form subproblem solutions

Step 4: Sorting and Clustering

Cell ordering: Hierarchical clustering on flattened motif profiles
Motif ordering: Hierarchical clustering on mean pairwise cross-correlation of occurrences

Both use optimal leaf ordering for visualization.

Improvements and Enhancements

From Original to Current MATLAB

The original implementation (918-line monolithic script) has been refactored into a proper MATLAB package (+tderica/):

Aspect	Original	Current
Structure	Single 918-line script with hardcoded paths	Modular package with 9 focused functions
Data Loading	Built-in Excel loading, normalization, outlier removal	Assumes pre-processed input
Path Handling	Hardcoded Windows path (`C:\Users\toyo\Desktop\...`)	Generic, user-provided input
Cache	Basic `.mat` file caching	Same, but cleaner naming conventions
Reproducibility	Fixed random seed behavior	Explicit `tderica_orders.mat` for paper figures
Documentation	Minimal inline comments	Full docstrings with input/output specs

Key functions extracted from the original:

doTdeRica -> tderica.decompose
doTdeRicaMatrixFactorization -> tderica.matrix_factorize
selectDataWithoutMissingValue -> tderica.select_subset
sortByClustering -> tderica.sort_cluster
makeFigures -> tderica.visualize
delayembed / delayembed_inv -> preserved verbatim
obtainCompEmbed -> tderica.project

From MATLAB to Python

The Python port adds significant algorithmic enhancements while maintaining numerical equivalence:

1. Parallel Time-Delay Embedding

from tderica import delayembed_batch

# Parallel embedding across all samples (joblib, auto-fallback to serial)
embed = delayembed_batch(data, dim_embed=300, n_jobs=-1)

2. Data-Driven Embedding Dimension (FNN)

Replaces the hard-coded dim_embed=300 with a principled estimate via the False Nearest Neighbours algorithm (Kennel et al., 1992):

from tderica import estimate_embed_dim

best_dim, fnn_fraction = estimate_embed_dim(signal, max_dim=500)
# best_dim: first dimension where false-neighbour fraction < 1%

3. ALS Solver for Matrix Factorization

Alternating Least Squares as an alternative to L-BFGS:

model = TdeRICA(dim_embed=300, n_components=14, method="als")
# or: matrix_factorize(..., method="als")

Trade-offs:

L-BFGS: Uses joint curvature information (faster for well-conditioned problems)
ALS: Guarantees monotone decrease via exact subproblem solves (more stable for ill-conditioned problems)

4. Cross-Validation (Leave-Cells-Out)

from tderica import cross_validate, select_n_components, select_dim_embed

cv = cross_validate(data, dim_embed=300, n_components=14, mask_ratio=0.2)
# cv["rmse"], cv["r2"], cv["correlation"]

best_n, results = select_n_components(data, dim_embed=300, n_candidates=[4, 8, 12, 14, 16, 20])
best_dim, results = select_dim_embed(data, n_components=14, dim_candidates=[100, 200, 300, 400])

5. Comprehensive Similarity Metrics Suite

Three-layer comparison framework for real vs simulated dynamics:

Layer 1 -- Time Alignment:

Function	Metric	Question Answered
`cosine_heatmap`	Cross-recurrence plot	When does pattern A appear in B?
`dtw_distance`	Dynamic Time Warping	Are trajectories similar (allowing warping)?
`frechet_distance`	Discrete Fréchet	What is the worst-case deviation?

Layer 2 -- Distribution:

Function	Metric	Question Answered
`wasserstein_distance`	Earth Mover's Distance	How different are state-space distributions?
`kl_divergence_1d`	Per-component KL	Which component's distribution mismatches?

Layer 3 -- Continuous Dynamics (D3):

Function	Method	Question Answered
`kernel_transition_comparison`	Gaussian-kernel Perron-Frobenius	Are transition probabilities consistent?
`transfer_entropy`	Kraskov-Stögbauer-Grassberger	How much does A's past inform B's future?
`local_jacobian_comparison`	k-NN local linearization	Are local expansion/rotation rates consistent?

6. Occurrence Analysis

from tderica import occurrence_cosine_similarity, sliding_window_similarity

sim = occurrence_cosine_similarity(occurrences, mode="between_samples")  # (num_comp, S, S)
sw  = sliding_window_similarity(occurrences, window_size=100, component=0)

7. Scikit-learn-style API

from tderica import TdeRICA

model = TdeRICA(dim_embed=300, n_components=14, lambda_=0.001)
model.fit(data)
occurrences, motifs = model.fit_transform(data)
rec = model.reconstruct()

Algorithmic Enhancements Summary

Enhancement	Original	Current MATLAB	Python
Parallel TDE	No	No	Yes (joblib)
FNN dim estimation	No	No	Yes
ALS solver	No	No	Yes
Cross-validation	No	No	Yes
Similarity metrics	Basic	Basic	Comprehensive (3-layer)
Occurrence analysis	No	No	Yes
scikit-learn API	No	No	Yes
FastICA fallback	No	No	Yes

Python Package Reference

Core Modules

Module	Key Functions	Description
`_tde.py`	`delayembed`, `delayembed_batch`, `delayembed_inv`, `estimate_embed_dim`	Time-delay embedding and dimension estimation
`_select_subset.py`	`select_subset`	Greedy NaN-free subset selection
`_project.py`	`project`	Project embedded data into motif space
`_rica.py`	`RICA`, `rica_fastica`	Reconstruction ICA solvers
`_decompose.py`	`decompose`	Full TDE-RICA pipeline
`_matrix_factorize.py`	`matrix_factorize`	L-BFGS or ALS imputation
`_sort_cluster.py`	`sort_clustering`	Hierarchical clustering and sorting
`_compare.py`	`compare_sim_to_real`	Real vs simulated comparison
`_validate.py`	`cross_validate`, `select_n_components`, `select_dim_embed`	Model selection
`_analysis.py`	`occurrence_cosine_similarity`, `sliding_window_similarity`, `temporal_profile_similarity`	Post-hoc analysis
`_similarity.py`	`cosine_heatmap`, `dtw_distance`, `wasserstein_distance`, `transfer_entropy`, `local_jacobian_comparison`, `similarity_report`	Comprehensive metrics
`_visualize.py`	`plot_sample_panel`, `plot_occurrence_scatter`, `plot_cross_correlation`, `plot_motifs`, `plot_pair_analysis`	Visualization

Similarity Metrics Suite

from tderica import similarity_report

report = similarity_report(comp_real, comp_sim, d3_fast=True)
# report["time_alignment"]   -- cosine stats, DTW, Frechet
# report["distribution"]    -- Wasserstein (global + per-component), KL
# report["dynamics"]        -- kernel transition, transfer entropy, Jacobian

MATLAB Package Reference

+tderica/
├── delayembed.m          -- Time-delay embedding
├── delayembed_inv.m      -- Inverse TDE (overlap-add average)
├── select_subset.m       -- Greedy NaN-free subset selection
├── decompose.m           -- TDE-RICA decomposition
├── matrix_factorize.m    -- L-BFGS matrix factorization for imputation
├── sort_cluster.m        -- Hierarchical clustering and sorting
├── visualize.m           -- Visualization functions
├── project.m             -- Project new data into motif space
└── compare_sim_to_real.m -- Compare real vs simulated trajectories

Biological Validation

Beyond numerical correctness, TDE-RICA outputs should be biologically reasonable. The following metrics validate the biological interpretability of extracted motifs and occurrences:

1. Motif Temporal Smoothness

Biological neural dynamics evolve smoothly in time. A motif with high roughness (large second temporal differences) suggests noise rather than a coherent dynamical pattern.

motif = motifs[k, :, :]  # (dim_embed, num_cells)
roughness = np.mean(np.var(np.diff(motif, n=2, axis=0), axis=1))
# Expect: roughness < 1.0 for biological data

2. Occurrence Temporal Autocorrelation

Motif activations should exhibit temporal structure. Near-zero or negative lag-1 autocorrelation suggests the occurrences are random or artifactual.

occ = occurrences[:, sample, k]  # (T_embed,)
acf = np.corrcoef(occ[1:], occ[:-1])[0, 1]
# Expect: acf > 0.3 for biological data

3. Motif Spatial Coherence

Cells participating in the same motif should have correlated spatial patterns. Anti-correlated or uncorrelated cells within a motif may indicate overfitting.

motif = motifs[k, :, :]  # (dim_embed, num_cells)
cormat = np.corrcoef(motif.T)
# Expect: mean pairwise correlation > 0.05

4. Component Orthogonality

ICA components should be approximately independent. Strong residual correlation between components suggests insufficient separation of dynamical modes.

comps = occurrences[:, sample, :]  # (T_embed, num_comp)
cormat = np.corrcoef(comps.T)
# Expect: mean |off-diagonal correlation| < 0.3

5. Reconstruction Fidelity

The model should capture a substantial fraction of data variance. Very low $R^2$ indicates underfitting or excessive noise.

# After decomposition
r2 = 1 - np.nansum((data - rec)**2) / np.nansum((data - np.nanmean(data))**2)
# Expect: R^2 > 0.5 for well-fit models

Validation with Real Data

When the .mat data file is available, run:

python -m pytest test_real_data_validation.py -v

This tests alignment with stored MATLAB results (16/20 animals within <1e-4 tolerance, all correlations >0.95) plus the five biological reasonableness metrics above.

Testing

Python Unit Tests

python3 tderica_test.py

14 tests covering TDE round-trip, subset selection, RICA, FastICA, projection, decomposition, matrix factorization, sorting, full pipeline, parallel TDE, FNN embedding dimension, ALS solver, cross-validation, and cosine similarity.

MATLAB-Python Alignment Tests

python -m pytest test_matlab_alignment.py -v

32 tests verifying numerical equivalence between MATLAB and Python implementations on controlled synthetic data, including:

TDE shape and value alignment
Inverse TDE round-trip and NaN handling
select_subset deterministic behavior
project least-squares consistency
RICA reconstruction quality
decompose shapes and reconstruction error
matrix_factorize full data extension
sort_clustering ordering consistency
cross-validate metric ranges
Edge cases (single sample, all NaN, empty embedding)
Numerical stability (large/small values)

Bugs Found and Fixed During Alignment

ALS solver stacking bug (_matrix_factorize.py): M_stack should use hstack (horizontal) not vstack (vertical) since M_base and M_add both have shape (num_comp, n_features). The target Y for step 1 should also use hstack([B, C]) not vstack.
Test data generation: The data_with_nan fixture originally used 15% random NaN which often made all cells invalid. Fixed by using structured NaN (different cells per sample) to ensure overlap.

Datasets

The original paper used whole-brain calcium imaging data from C. elegans:

Citation

If you use this code, please cite:

Toyoshima, Y., Sato, H., Nagata, D., Kanamori, M., Jang, M. S., Kuze, K., Oe, S., Teramoto, T., Iwasaki, Y., Ishihara, T., & Iino, Y. (2024). Ensemble dynamics and information flow deduction from whole-brain imaging data. PLOS Computational Biology, 20(3), e1011848. https://doi.org/10.1371/journal.pcbi.1011848

License

MIT License. See LICENSE for details.

Acknowledgments

Original MATLAB implementation: Yuichi Toyoshima (University of Tokyo)
Python port and enhancements: Project MAGI
The False Nearest Neighbours implementation follows Kennel et al. (1992)
The RICA objective follows the MATLAB Statistics and Machine Learning Toolbox specification

Name		Name	Last commit message	Last commit date
+tderica		+tderica
octave_compat		octave_compat
tderica		tderica
.gitignore		.gitignore
BIOLOGY_MATH_ANALYSIS.md		BIOLOGY_MATH_ANALYSIS.md
LICENSE		LICENSE
README.md		README.md
demo_compare_real_and_simulation.m		demo_compare_real_and_simulation.m
setup.py		setup.py
setup_octave.m		setup_octave.m
tderica_main.m		tderica_main.m
tderica_orders.mat		tderica_orders.mat
tderica_test.py		tderica_test.py
test_matlab_alignment.py		test_matlab_alignment.py
test_real_data_comparison.m		test_real_data_comparison.m
test_real_data_validation.py		test_real_data_validation.py

Latest commit

History

Folders and files

Repository files navigation