diff --git a/README.md b/README.md index 731e8fb..49af894 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,8 @@ GPLA summarizes high-dimensional **spike–[Local Field Potential (LFP)](https:/ - **Statistical testing**: - fast **RMT-based** heuristic (Marchenko–Pastur edge) - **spike-jitter** surrogate tests (interval / ISI-preserved / group-preserved / population) -- **Preprocessing hooks**: trial concatenation, spike-count filtering, optional PCA whitening +- **Data preparation**: trial concatenation, spike-count filtering, temporal/unit selection, + and optional PCA whitening of user-provided analytic LFP signals - **Simulations**: phase-locked and transient coupling generators - **Figure reproduction**: a Figure 2-style simulation and visualization script @@ -106,7 +107,13 @@ print("Spike vector shape:", result.spike_vector.shape) Input conventions: - spike trains: list of trials, each `(n_units, n_samples)` -- analytic LFP: `(n_channels, n_samples, n_trials)` complex array (bandpass + Hilbert in real data) +- preferred LFP input: `(n_channels, n_samples, n_trials)` complex analytic array +- supported alternative: a real array of phase angles in radians + +> `gpla()` does not treat real-valued input as raw LFP voltage. Raw LFP data must be +> band-pass filtered and converted to an analytic signal upstream, for example with a +> Hilbert transform. If a real array is passed, PyGPLA warns that it will be interpreted +> as phase angles in radians. ## Reproducing Figure 2 from the original paper (Safavi et al., 2023) @@ -141,17 +148,45 @@ pytest ## Citing -If you use PyGPLA in your research, please cite the original GPLA method paper: +If you use PyGPLA in your research, please cite the following papers. + +The original GPLA method paper: -> Safavi, S., et al. (2023). *Uncovering the organization of neural circuits with -> Generalized Phase Locking Analysis.* PLOS Computational Biology. +> Safavi, S., Panagiotaropoulos, T. I., Kapoor, V., Ramirez-Villegas, J. F., +> Logothetis, N. K., & Besserve, M. (2023). *Uncovering the organization of neural +> circuits with Generalized Phase Locking Analysis.* PLOS Computational Biology, +> 19(4), e1010983. ```bibtex @article{safavi2023uncovering, title = {Uncovering the organization of neural circuits with Generalized Phase Locking Analysis}, - author = {Safavi, Shervin and others}, - journal = {PLOS Computational Biology}, - year = {2023} + author = {Safavi, Shervin and Panagiotaropoulos, Theofanis I. and Kapoor, Vishal and Ramirez-Villegas, Juan F. and Logothetis, Nikos K. and Besserve, Michel}, + journal = {PLoS Computational Biology}, + year = {2023}, + volume = {19}, + number = {4}, + pages = {e1010983}, + doi = {10.1371/journal.pcbi.1010983} +} +``` + +The mathematical framework underlying GPLA's coupling measure and its Random Matrix +Theory–based significance test: + +> Safavi, S., Logothetis, N. K., & Besserve, M. (2021). *From Univariate to +> Multivariate Coupling Between Continuous Signals and Point Processes: A +> Mathematical Framework.* Neural Computation, 33(7), 1751–1817. + +```bibtex +@article{safavi2021univariate, + title = {From Univariate to Multivariate Coupling Between Continuous Signals and Point Processes: A Mathematical Framework}, + author = {Safavi, Shervin and Logothetis, Nikos K. and Besserve, Michel}, + journal = {Neural Computation}, + year = {2021}, + volume = {33}, + number = {7}, + pages = {1751--1817}, + doi = {10.1162/neco_a_01389} } ``` diff --git a/docs/quickstart.md b/docs/quickstart.md index b55a122..69fbd39 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -116,12 +116,14 @@ print("Selected units:", result.metadata["selected_units"]) `gpla(spike_trains, lfp_signal, ...)` expects: - `spike_trains`: list of length `n_trials`, each element a 2D array `(n_units, n_samples)` -- `lfp_signal`: a 3D **complex** array `(n_channels, n_samples, n_trials)` +- `lfp_signal`: preferably a 3D **complex analytic** array + `(n_channels, n_samples, n_trials)` -:::{note} -For real data, you typically build `lfp_signal` by bandpass filtering around a target -frequency and applying a Hilbert transform to get the complex analytic signal. -The Figure 2 tutorial shows a concrete example of that preprocessing: {doc}`tutorials`. +:::{important} +`gpla()` does not interpret a real-valued array as raw LFP voltage. Real input is supported +only as phase angles in radians and produces a warning. For raw LFP data, first bandpass +filter around the frequency band of interest and apply a Hilbert transform to obtain the +complex analytic signal. The Figure 2 tutorial shows a concrete example: {doc}`tutorials`. ::: ## Interpreting the outputs (what you get back) @@ -171,6 +173,10 @@ def bandpass_hilbert(x: np.ndarray, sf: float, band_hz: tuple[float, float]) -> You can then pass `lfp_analytic` into `{py:func}`pygpla.api.gpla``. +If you already have phase angles in radians, you may pass that real-valued phase array +directly. PyGPLA emits a warning to distinguish this supported phase representation from +accidentally supplied raw LFP voltage. + ## Adding significance testing (optional) PyGPLA supports two main significance-testing modes: diff --git a/docs/software_architecture.md b/docs/software_architecture.md index 28509a4..6984e11 100644 --- a/docs/software_architecture.md +++ b/docs/software_architecture.md @@ -110,9 +110,12 @@ This keeps the numerical path explicit and lightweight: The package assumes: - spikes as trial list of `(units, samples)` arrays -- analytic LFP as `(channels, samples, trials)` complex array +- preferably an analytic LFP as `(channels, samples, trials)` complex array +- alternatively, phase angles in radians as a real array of the same shape -These contracts are enforced in preprocessing and validation utilities. +Raw LFP voltage must be converted to the desired frequency-specific analytic signal upstream. +The public API warns when it receives real input so that raw voltage is not silently confused +with the supported phase-angle representation. ### Deterministic and Stochastic Paths diff --git a/docs/usage.md b/docs/usage.md index 0271821..786b15d 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -24,7 +24,10 @@ pip install -e .[docs] ## Data expectations - **Spike trains**: list of arrays shaped `(units, samples)` per trial, binary or counts. -- **LFP analytic signal**: complex array shaped `(channels, samples, trials)`; use a bandpass + Hilbert transform upstream. +- **LFP analytic signal (preferred)**: complex array shaped `(channels, samples, trials)`; + use a bandpass filter and Hilbert transform upstream when starting from raw LFP voltage. +- **LFP phase representation (supported)**: real array of the same shape containing phase + angles in radians. Real input triggers a warning because it is not interpreted as raw voltage. - **Sampling**: spike and LFP sample counts must match within each trial; supply `sampling_frequency` when using jitter surrogates. ## Minimal example @@ -65,6 +68,9 @@ print("Spike vector shape:", result.spike_vector.shape) ## Preprocessing knobs +These options prepare an analytic-signal or phase input for GPLA; they do not bandpass-filter +raw LFP voltage or compute its Hilbert transform. + - `plvNrmlzMethed`: `nSpk`, `nSpk-square-root` (default), or `var1_theoretical`. - `flag_whitening`: 0 (off), 1/2 for PCA whitening variants; optionally set `PreprocessingConfig.whitening.variance_proportion`. - `flag_lfpNrmlz`: normalize analytic LFP amplitude if set. diff --git a/paper/paper.md b/paper/paper.md index f40d72d..41200e1 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -37,7 +37,7 @@ url: "https://pygpla.readthedocs.io/en/latest/" ## Summary -PyGPLA is a Python implementation of Generalized Phase Locking Analysis (GPLA) for multivariate analysis of coupling between spikes and local field potentials (LFPs) [@safavi2023uncovering]. For a given frequency, GPLA constructs a complex coupling matrix $\hat{C}(f) \in \mathbb{C}^{N_c \times N_u}$ between LFP channels ($N_c$) and spike units ($N_u$), then applies singular value decomposition (SVD) to reduce the dimensionality of data. The leading singular value summarizes population-level coupling strength, while the corresponding singular vectors describe dominant LFP and spike coupling modes. PyGPLA provides a complete analysis workflow including LFP preprocessing, coupling-matrix construction, SVD-based decomposition, and diverse statistical significance testing [@safavi2021univariate]. +PyGPLA is a Python implementation of Generalized Phase Locking Analysis (GPLA) for multivariate analysis of coupling between spikes and local field potentials (LFPs) [@safavi2023uncovering]. For a given frequency, GPLA constructs a complex coupling matrix $\hat{C}(f) \in \mathbb{C}^{N_c \times N_u}$ between LFP channels ($N_c$) and spike units ($N_u$), then applies singular value decomposition (SVD) to reduce the dimensionality of data. The leading singular value summarizes population-level coupling strength, while the corresponding singular vectors describe dominant LFP and spike coupling modes. PyGPLA accepts a user-provided frequency-specific analytic LFP signal or phase representation and provides data selection, optional whitening and normalization, coupling-matrix construction, SVD-based decomposition, and statistical significance testing [@safavi2021univariate]. ## Statement of need @@ -56,7 +56,7 @@ Existing spike–field coupling methods - including the phase-locking value (PLV PyGPLA provides a comprehensive solution for multivariate spike–field coupling analysis, including: -- **LFP preprocessing:** Band-pass filtering, Hilbert transform for analytic signal extraction [@chavez2006proper], and optional reduced-rank whitening to decorrelate channels while avoiding noise amplification [@safavi2023uncovering]. +- **LFP input and data preparation:** PyGPLA accepts a user-provided complex analytic LFP signal or a real phase representation in radians. Raw LFP voltage must first be converted upstream to the desired frequency-specific analytic signal, for example through band-pass filtering followed by a Hilbert transform [@chavez2006proper]. PyGPLA then supports temporal and unit selection, optional channel-wise normalization, and optional reduced-rank whitening to decorrelate channels while avoiding noise amplification [@safavi2023uncovering]. - **Coupling-matrix construction:** Assembly of the complex-valued coupling matrix $\widehat{\mathbf{C}}(f) \in \mathbb{C}^{N_c \times N_u}$, where each entry sums the analytic LFP evaluated at all spike times of a given unit. - **SVD-based decomposition:** Extraction of the generalized phase locking value or gPLV (the leading singular value) and associated LFP and spike spatial vectors, with rotational phase alignment, unwhitening of LFP vectors, and spike-vector rescaling. - **Statistical testing:** Significance assessment via two complementary approaches: (1) surrogate-based testing using multiple spike-jittering schemes, and (2) an analytical test based on Marchenko–Pastur Random Matrix Theory (RMT) [@anderson2010random; @safavi2023uncovering]. diff --git a/src/pygpla/api.py b/src/pygpla/api.py index f93b35c..1ce801d 100644 --- a/src/pygpla/api.py +++ b/src/pygpla/api.py @@ -2,6 +2,7 @@ from __future__ import annotations +import warnings from dataclasses import asdict, dataclass, is_dataclass from typing import Any, Dict, Optional, Union @@ -89,7 +90,9 @@ def gpla( spike_trains : List of spike arrays shaped (units, samples) per trial. lfp_signal : - Complex analytic LFP array shaped (channels, samples, trials). + Complex analytic LFP array shaped (channels, samples, trials). Real-valued + inputs are interpreted as phase angles in radians, not as raw LFP voltage, + and produce a warning. flag_gPLVnrmlz : Legacy gPLV normalization flag (0 keep raw, nonzero scales by matrix size). nSpikeThreshold : @@ -140,6 +143,16 @@ def gpla( spike_list = [np.asarray(st) for st in spike_trains] lfp_array = np.asarray(lfp_signal) + if np.isrealobj(lfp_array): + warnings.warn( + "gpla received a real-valued LFP array. Real inputs are interpreted " + "as phase angles in radians, not as raw LFP voltage. If this is raw " + "LFP data, band-pass filter it and compute its analytic signal " + "(for example, using a Hilbert transform) before calling gpla().", + UserWarning, + stacklevel=2, + ) + ( spikeTrains_allTrLong, lfpPhases_allTrLong, diff --git a/src/pygpla/config.py b/src/pygpla/config.py index 009c8a2..3269c83 100644 --- a/src/pygpla/config.py +++ b/src/pygpla/config.py @@ -69,10 +69,11 @@ def validate_spike_trains(spike_trains: Sequence[np.ndarray]) -> None: def validate_lfp_signal(lfp_signal: ArrayLike) -> np.ndarray: """ - Ensure LFP input can be treated as a complex analytic array. + Validate the dimensions of an analytic-LFP or phase-array input. Requires a 3D array shaped (channels, samples, trials); raises ValueError otherwise. - Returns the array as np.ndarray for downstream use. + Complex arrays represent analytic LFP signals. Real arrays are interpreted downstream + as phase angles in radians, not as raw LFP voltage. Returns the input as an ndarray. """ arr = np.asarray(lfp_signal) diff --git a/src/pygpla/core/coupling.py b/src/pygpla/core/coupling.py index 60b00be..fa5f03b 100644 --- a/src/pygpla/core/coupling.py +++ b/src/pygpla/core/coupling.py @@ -56,8 +56,10 @@ def compute_coupling_matrix( spike_trains: Array of shape (units, samples) with binary or count spikes. lfp_signal: - Real or complex LFP array of shape (channels, samples). Real inputs are wrapped - to unit-magnitude complex exponentials; complex inputs preserve amplitude. + Real phase-angle or complex analytic-LFP array of shape (channels, samples). + Real inputs are interpreted as radians and wrapped to unit-magnitude complex + exponentials; they are not treated as raw LFP voltage. Complex inputs preserve + amplitude. normalization_method: One of ``{"nSpk", "nSpk-square-root", "var1_theoretical"}`` controlling spike-count scaling in the coupling matrix. diff --git a/tests/test_api.py b/tests/test_api.py new file mode 100644 index 0000000..9eea4ca --- /dev/null +++ b/tests/test_api.py @@ -0,0 +1,35 @@ +"""Tests for the public GPLA API.""" + +import warnings + +import numpy as np +import pytest + +from pygpla.api import gpla + + +def _minimal_inputs(): + spikes = [np.array([[0, 1, 0, 0]], dtype=int)] + lfp_analytic = np.ones((1, 4, 1), dtype=complex) + return spikes, lfp_analytic + + +def test_gpla_warns_that_real_lfp_is_interpreted_as_phase(): + spikes, lfp_analytic = _minimal_inputs() + lfp_phase = np.angle(lfp_analytic) + + with pytest.warns(UserWarning, match="phase angles in radians, not as raw LFP voltage"): + result = gpla(spikes, lfp_phase, plvNrmlzMethed="nSpk") + + assert np.isfinite(result.gplv) + + +def test_gpla_does_not_warn_for_complex_analytic_lfp(): + spikes, lfp_analytic = _minimal_inputs() + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + result = gpla(spikes, lfp_analytic, plvNrmlzMethed="nSpk") + + assert not caught + assert np.isfinite(result.gplv)