diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml index 01d84a7..7ea4c72 100644 --- a/.github/workflows/draft-pdf.yml +++ b/.github/workflows/draft-pdf.yml @@ -11,6 +11,9 @@ on: - ".github/workflows/draft-pdf.yml" workflow_dispatch: +permissions: + contents: read + jobs: paper: runs-on: ubuntu-latest @@ -33,3 +36,4 @@ jobs: # This is the output path where Pandoc will write the compiled # PDF. Note, this should be the same directory as the input paper.md path: paper/paper.pdf + if-no-files-found: error diff --git a/paper/paper.md b/paper/paper.md index c12f393..65d88bf 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -11,23 +11,31 @@ authors: - name: Amir Khani orcid: 0009-0008-8042-9590 corresponding: true - affiliation: 1 + affiliation: "1, 2" + - name: Christian Beste + affiliation: "3, 4" - name: Shervin Safavi orcid: 0000-0002-2868-530X corresponding: true - affiliation: "2, 3" + affiliation: "5, 6" affiliations: - name: Donders Institute for Brain, Cognition, and Behaviour, Radboud University, Nijmegen, Netherlands index: 1 ror: "016xsfp80" - - name: Computational Neuroscience, Department of Child and Adolescent Psychiatry, Faculty of Medicine, Technische Universität Dresden, Dresden 01307, Germany + - name: Department of Computer Science, Amirkabir University of Technology, Tehran, Iran index: 2 + - name: Cognitive Neurophysiology, Department of Child and Adolescent Psychiatry, Faculty of Medicine, Technische Universität Dresden, Dresden, Germany + index: 3 + - name: German Center for Child and Adolescent Health (DZKJ), partner site Leipzig/Dresden, Dresden, Germany + index: 4 + - name: Computational Neuroscience, Department of Child and Adolescent Psychiatry, Faculty of Medicine, Technische Universität Dresden, Dresden 01307, Germany + index: 5 ror: "042aqky30" - name: Department of Computational Neuroscience, Max Planck Institute for Biological Cybernetics, Tübingen 72076, Germany - index: 3 + index: 6 ror: "026nmvv73" -date: "2026-03-22" +date: "22 March 2026" bibliography: paper.bib repository: "https://github.com/CMC-lab/PyGPLA" crossref: true @@ -48,7 +56,12 @@ However, commonly used spike-LFP coupling measures are pairwise, making them sub GPLA was developed to address these challenges by providing an efficient multivariate framework together with statistical routines [@safavi2021univariate] for characterizing spike-LFP coupling at the population level [@safavi2023uncovering]. However, the original algorithm was implemented in MATLAB, limiting its accessibility. PyGPLA bridges this gap as an open-source Python implementation, aligning with the extensive use of Python in neuroscience as solidified by libraries such as MNE-Python [@GramfortEtAl2013a], Nilearn [@Nilearn], and TranCIT [@Nouri2025]. -Existing spike–field coupling methods - including the phase-locking value (PLV), pairwise phase consistency (PPC) [@vinck2010pairwise], and spike-field coherence - operate on individual spike–LFP channel pairs. While informative for small-scale recordings, these pairwise approaches do not directly capture the population-level structure of coupling, which is of paramount importance for recordings with modern high-density probes. GPLA addresses this gap by providing a multivariate method that summarizes all spike–field interactions simultaneously, analogous to how multivariate statitical methods summarizes covariance structure (well-known example is Principal component analysis or PCA). To our knowledge, no other Python package implements this multivariate approach to spike–field coupling analysis. +Existing spike–field coupling methods - including the phase-locking value (PLV), pairwise phase consistency (PPC) [@vinck2010pairwise], and spike-field coherence - operate on individual spike–LFP channel pairs. While informative for small-scale recordings, these pairwise approaches do not directly capture the population-level structure of coupling, which is of paramount importance for recordings with modern high-density probes. GPLA addresses this gap by providing a multivariate method that summarizes all spike–field interactions simultaneously, analogous to how multivariate statistical methods such as principal component analysis (PCA) summarize covariance structure. To our knowledge, no other Python package implements this multivariate approach to spike–field coupling analysis. + +## State of the field + +Common spike–field coupling measures, including the phase-locking value, pairwise phase consistency, and spike–field coherence, quantify relationships between individual spike units and LFP channels [@vinck2010pairwise; @vinck2012improved]. General-purpose Python packages such as MNE-Python provide extensive signal-processing and electrophysiology infrastructure [@GramfortEtAl2013a], but do not implement GPLA’s joint decomposition of the complete spike–field coupling matrix. The only directly comparable implementation is the MATLAB research code accompanying the original GPLA study [@safavi2023uncovering]. +PyGPLA was developed as a separate package because GPLA requires a method-specific pipeline comprising coupling-matrix construction, normalization, optional reduced-rank whitening, SVD-based extraction of population coupling modes, and analytical or surrogate-based statistical inference. These operations do not extend an existing pairwise measure or general preprocessing function in a natural way. PyGPLA’s scholarly contribution is therefore a dedicated, open-source Python implementation of the complete GPLA methodology. ## Functionality @@ -64,24 +77,33 @@ PyGPLA provides a comprehensive solution for multivariate spike–field coupling ### Normalization options -The package supports two normalization modes for the coupling matrix: PLV-type normalization ($1/N_m^{\mathrm{tot}}$) for phase-only analysis, and unit-variance normalization ($1/\sqrt{N_m^{\mathrm{tot}}}$) required for the theoretical RMT-based significance test. +The package supports two normalization modes for the coupling matrix, expressed in terms of the total spike count $N_m^{\mathrm{tot}}$ for unit $m$: PLV-type normalization ($1/N_m^{\mathrm{tot}}$) for phase-only analysis, and unit-variance normalization ($1/\sqrt{N_m^{\mathrm{tot}}}$) required for the theoretical RMT-based significance test. ## Example We illustrate PyGPLA on synthetic transient-coupling simulations generated with the script `paper/figures/figure2.py`. As shown in \autoref{fig:gpla_results}, the figure compares four coupling models and reports their recovered gPLV values, together with representative LFP and spike-train windows and the associated spike-vector structure. -![Illustrative PyGPLA simulation. (A) LFP analysis window. (B) gPLV across models M1-M4. (C-F) Model schematics, LFP and spikes, and recovered spike vectors.](figures/figure2_python.png){#fig:gpla_results width="70%"} +![Illustrative PyGPLA simulation. (A) LFP analysis window. (B) gPLV across Models 1-4. (C-F) Model schematics, LFP and spikes, and recovered spike vectors.](figures/figure2_python.png){#fig:gpla_results width="70%"} Code snippets and detailed instructions for reproducing these results are available in the package documentation and example scripts in the repository. -## Implementation details +## Software design + +PyGPLA uses a function-oriented, layered design that separates the main stages of the analysis while providing a unified high-level workflow. The high-level `gpla()` function coordinates data preparation, GPLA decomposition, and optional statistical testing, returning the results and relevant bookkeeping in a single `GPLAResult` object. The underlying operations—coupling-matrix construction, SVD factorization, whitening, jitter generation, and simulation—remain independently accessible. This design provides a concise default workflow while allowing researchers to inspect, test, or replace individual methodological stages. Automated tests validate these independently accessible numerical components. + +PyGPLA accepts standard NumPy arrays rather than requiring a package-specific data container, facilitating integration with existing electrophysiology workflows. Frequency selection and conversion of raw LFP voltage to an analytic signal are intentionally left upstream because these operations require experiment-specific filtering choices. PyGPLA therefore operates on a frequency-specific complex analytic signal or phase representation and warns when real-valued input is supplied. + +Several numerical and interface conventions preserve continuity with the original MATLAB implementation, facilitating validation against the reference implementation and migration of existing GPLA analyses. These include its normalization alternatives, phase convention, reduced-rank whitening methods, and selected legacy parameter names. Optional PCA-based whitening reduces correlations among LFP channels, while an unwhitening operator maps the resulting coupling modes back to the original channel coordinates. Statistical inference is separated from the deterministic decomposition, allowing users to choose between a computationally inexpensive analytical RMT-based decision and more expensive spike-jitter surrogate tests. The core package depends only on NumPy; SciPy is provided as an optional dependency for documented simulation and signal-processing workflows. + +## Research impact statement -The `pygpla` package is distributed under the BSD-2-Clause license. PyGPLA follows a layered architecture separating preprocessing (`pygpla.preprocessing`), core computation (`pygpla.core`), statistical testing (`pygpla.stats`), simulation utilities (`pygpla.simulations`), and configuration (`pygpla.config`). The high-level API (`pygpla.api.gpla`) orchestrates the full pipeline in a single call, keeping default usage compact while preserving access to lower-level components for advanced analyses. The package is built on NumPy [@harris2020array] and SciPy [@virtanen2020fundamental], and includes a `pytest` test suite with continuous integration via GitHub Actions. +GPLA has been evaluated in published simulations, biophysical network models, and multielectrode recordings, where it revealed population-level spike–field coupling patterns related to properties of the underlying neural circuits [@safavi2023uncovering]. PyGPLA transfers this established methodology from MATLAB research code into an installable, open-source Python package. To our knowledge, it is the first Python implementation that jointly decomposes the complete spike–field coupling matrix. Its immediate research contribution is therefore to make population-level GPLA available to researchers using Python-based electrophysiology workflows. Because PyGPLA is newly released, broader external adoption and independent applications remain to be established. -## Use of generative AI +## AI usage disclosure -Generative AI tools were used only to assist with drafting documentation and -refactoring code for style. All AI-assisted output was reviewed, tested, and +Generative artificial intelligence (AI) tools were used to assist with drafting +documentation and portions of this paper, and with refactoring code for style. +All AI-assisted output was reviewed, tested, and validated by the authors, who take full responsibility for the correctness of the software and the content of this paper. The scientific method, algorithmic design, and numerical implementation of GPLA were carried out by the authors. diff --git a/paper/paper.pdf b/paper/paper.pdf index b9bb7ab..c174885 100644 Binary files a/paper/paper.pdf and b/paper/paper.pdf differ diff --git a/src/pygpla/api.py b/src/pygpla/api.py index 1ce801d..24052af 100644 --- a/src/pygpla/api.py +++ b/src/pygpla/api.py @@ -3,7 +3,7 @@ from __future__ import annotations import warnings -from dataclasses import asdict, dataclass, is_dataclass +from dataclasses import dataclass from typing import Any, Dict, Optional, Union import numpy as np @@ -59,8 +59,19 @@ def _stats_config_to_dict( return None if isinstance(stats_config, dict): return stats_config - if is_dataclass(stats_config): - return asdict(stats_config) + if isinstance(stats_config, StatTestConfig): + legacy_config = dict(stats_config.additional_params) + legacy_config.update( + { + "testType": stats_config.test_type, + "nJtr": stats_config.n_surrogates, + "jitterWinWidth": stats_config.jitter_window_width, + "spkSF": stats_config.sampling_frequency, + "jitterType": stats_config.jitter_type, + "alphaValue": stats_config.alpha, + } + ) + return legacy_config raise TypeError("stats_config must be a dict, StatTestConfig, or None.") @@ -153,6 +164,8 @@ def gpla( stacklevel=2, ) + stat_dict = _stats_config_to_dict(stats_config) + ( spikeTrains_allTrLong, lfpPhases_allTrLong, @@ -167,7 +180,7 @@ def gpla( unitSubset=unitSubset, temporalWindow=temporalWindow, flag_origDimEigVec=flag_origDimEigVec, - statTestInfo=stats_config, + statTestInfo=stat_dict, iSV=iSV, checkSameElecStuff_flag=0, plvNrmlzMethed=plvNrmlzMethed, @@ -187,8 +200,6 @@ def gpla( else: sameElecCheckInfo = None - stat_dict = _stats_config_to_dict(stats_config) - ( gPLV, pValue, diff --git a/src/pygpla/config.py b/src/pygpla/config.py index 3269c83..10657ef 100644 --- a/src/pygpla/config.py +++ b/src/pygpla/config.py @@ -3,7 +3,7 @@ from __future__ import annotations from dataclasses import dataclass, field -from typing import Dict, Optional, Sequence, Tuple, Union +from typing import Any, Dict, Optional, Sequence, Tuple, Union import numpy as np @@ -37,7 +37,12 @@ class PreprocessingConfig: @dataclass(slots=True) class StatTestConfig: - """Statistical testing controls for the GPLA pipeline.""" + """Statistical testing controls for the GPLA pipeline. + + ``additional_params`` accepts specialized legacy runner options not represented + by named fields, such as ``rngSeed`` or ``SVspectrumStatsType``. Named fields + take precedence if the same legacy key is also present in ``additional_params``. + """ test_type: str = "RMT-based" n_surrogates: int = 0 @@ -45,7 +50,7 @@ class StatTestConfig: sampling_frequency: float = 1.0 jitter_type: str = "interval-jittering" alpha: float = 0.05 - additional_params: Dict[str, ArrayLike] = field(default_factory=dict) + additional_params: Dict[str, Any] = field(default_factory=dict) def validate_spike_trains(spike_trains: Sequence[np.ndarray]) -> None: diff --git a/tests/test_api.py b/tests/test_api.py index 9eea4ca..1be7bbc 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -5,7 +5,8 @@ import numpy as np import pytest -from pygpla.api import gpla +from pygpla.api import _stats_config_to_dict, gpla +from pygpla.config import StatTestConfig def _minimal_inputs(): @@ -33,3 +34,67 @@ def test_gpla_does_not_warn_for_complex_analytic_lfp(): assert not caught assert np.isfinite(result.gplv) + + +def test_gpla_accepts_stat_test_config_for_rmt(): + spikes, lfp_analytic = _minimal_inputs() + + result = gpla( + spikes, + lfp_analytic, + stats_config=StatTestConfig(test_type="RMT-based"), + plvNrmlzMethed="var1_theoretical", + ) + + assert np.isnan(result.p_value) + assert isinstance(result.stats["gPLV_stats"]["nullHypoReject"], bool) + + +def test_gpla_accepts_stat_test_config_for_spike_jittering(): + spikes, lfp_analytic = _minimal_inputs() + config = StatTestConfig( + test_type="spike-jittering", + n_surrogates=2, + jitter_window_width=0.01, + sampling_frequency=1000.0, + jitter_type="fake-jittering", + alpha=0.01, + additional_params={"rngSeed": 123, "SVspectrumStatsType": "default"}, + ) + + result = gpla( + spikes, + lfp_analytic, + stats_config=config, + plvNrmlzMethed="nSpk", + ) + + assert result.p_value == pytest.approx(0.0) + assert result.stats["gPLV_stats"]["nullDistribution"].shape == (2,) + + +def test_stat_test_config_maps_all_fields_to_legacy_keys(): + config = StatTestConfig( + test_type="spike-jittering", + n_surrogates=25, + jitter_window_width=0.02, + sampling_frequency=500.0, + jitter_type="group-preserved-interval-jittering", + alpha=0.01, + additional_params={ + "rngSeed": 42, + "SVspectrumStatsType": "RMT-heuristic", + "testType": "ignored-in-favor-of-named-field", + }, + ) + + assert _stats_config_to_dict(config) == { + "testType": "spike-jittering", + "nJtr": 25, + "jitterWinWidth": 0.02, + "spkSF": 500.0, + "jitterType": "group-preserved-interval-jittering", + "alphaValue": 0.01, + "rngSeed": 42, + "SVspectrumStatsType": "RMT-heuristic", + }