Joss review fixes

.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

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.

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,

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,15 +37,20 @@ 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
     jitter_window_width: float = 0.05
     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:

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",
+    }