Phase-Amplitude Coupling — Cross-Frequency Interaction¶

The pac module measures phase-amplitude coupling (PAC): the statistical dependence between the phase of a low-frequency oscillation and the amplitude (envelope) of a high-frequency oscillation. PAC is a hallmark of cross-frequency communication in neural circuits and a key biomarker in neuroscience.

In the SCPN framework, PAC quantifies how one layer's phase dynamics modulate another layer's amplitude — capturing hierarchical information flow across timescales.

1. Mathematical Formalism¶

1.1 The Modulation Index (MI)¶

The modulation index (Tort et al. 2010) quantifies how non-uniformly high-frequency amplitude is distributed across low-frequency phase bins:

Divide the low-frequency phase \(\theta_{\text{low}}\) into \(B\) equal bins on \([0, 2\pi)\): \([\phi_0, \phi_1), [\phi_1, \phi_2), \ldots\)
For each bin \(k\), compute the mean high-frequency amplitude: \(\bar{A}_k = \langle A_{\text{high}} \rangle_{\theta \in [\phi_k, \phi_{k+1})}\)
Normalise to a probability distribution: \(p_k = \bar{A}_k / \sum_k \bar{A}_k\)
Compute KL divergence from uniform: \(D_{KL} = \sum_k p_k \log(p_k \cdot B)\)
Normalise: \(\text{MI} = D_{KL} / \log(B)\)

\[ \text{MI} = \frac{1}{\log B} \sum_{k=1}^{B} p_k \log(B \cdot p_k) \in [0, 1] \]

MI value	Interpretation
MI = 0	No coupling — amplitude uniform across phase bins
0 < MI < 0.3	Weak coupling
0.3 ≤ MI < 0.6	Moderate coupling
MI ≥ 0.6	Strong coupling — amplitude strongly modulated by phase
MI = 1	Perfect coupling — all amplitude concentrated in one bin

1.2 Information-Theoretic Interpretation¶

The MI is a normalised KL divergence from the uniform distribution \(q_k = 1/B\):

\[ \text{MI} = \frac{D_{KL}(p \| q)}{\log B} = \frac{H_{\max} - H(p)}{H_{\max}} \]

where \(H(p) = -\sum p_k \log p_k\) is the entropy of the amplitude distribution across phase bins and \(H_{\max} = \log B\) is the entropy of the uniform distribution.

MI = 0 when entropy is maximal (uniform) and MI = 1 when entropy is zero (delta function — all amplitude in one bin).

1.3 PAC Matrix¶

For \(N\) oscillators with both phase \(\theta_i(t)\) and amplitude \(A_i(t)\) time series, the PAC matrix is:

\[ \text{PAC}_{ij} = \text{MI}(\theta_i, A_j) \]

Entry \((i, j)\) measures how much oscillator \(i\)'s phase modulates oscillator \(j\)'s amplitude. The matrix is generally NOT symmetric — the low-frequency phase driver and high-frequency amplitude modulated signal play different roles.

Diagonal entries \(\text{PAC}_{ii}\) measure self-coupling (within-signal PAC). Off-diagonal entries \(\text{PAC}_{ij}, i \neq j\) measure cross-signal PAC (inter-layer communication in SCPN).

1.4 PAC Gate¶

The pac_gate(pac_value, threshold) function provides a binary decision: is PAC significant? Default threshold 0.3 is based on empirical studies showing that MI > 0.3 reliably indicates genuine coupling above noise floor for typical neuroscience data (Tort et al. 2010).

1.5 Number of Phase Bins¶

The choice of \(B\) (default 18 = 20° bins) affects MI estimation:

B	Bin width	Resolution	Bias	Variance
6	60°	Coarse	Low	Low
18	20°	Standard	Medium	Medium
36	10°	Fine	High	High

Bias-variance trade-off: Fewer bins → lower variance but may miss narrow coupling peaks. More bins → higher resolution but increased noise (many bins may have few samples). Rule of thumb: \(B = 18\) for \(T > 500\) samples; reduce to \(B = 9\) for \(T < 200\).

1.6 Statistical Significance¶

MI is biased upward for finite data — even random signals produce MI > 0. To assess significance:

Surrogate testing: Shuffle amplitude time series, recompute MI. Repeat 200 times. Genuine MI should exceed 95th percentile of surrogates.
Analytical bound: For \(T\) i.i.d. samples and \(B\) bins: \(\text{MI}_{\text{null}} \sim \frac{B-1}{2T \ln B}\)
Permutation z-score: \(z = (\text{MI} - \mu_{\text{surr}}) / \sigma_{\text{surr}}\). Significant if \(z > 3\).

The module does NOT perform surrogate testing automatically. Users should implement it for publication-quality results.

2. Theoretical Context¶

2.1 Historical Background¶

Phase-amplitude coupling was first systematically characterised in neuroscience by Canolty et al. (2006) who observed theta (4–8 Hz) phase modulating gamma (30–100 Hz) amplitude in human electrocorticography. The modulation index method was formalised by Tort, Komorowski, Eichenbaum & Kopell (2010), becoming the standard PAC measure.

Jensen & Colgin (2007) proposed that PAC is the mechanism by which slow rhythms coordinate fast local processing — the "communication through coherence" framework.

2.2 Neuroscience Applications¶

Brain region	Low freq	High freq	Function
Hippocampus	Theta (4–8 Hz)	Gamma (30–100 Hz)	Memory encoding
Cortex	Alpha (8–12 Hz)	High-gamma (60–150 Hz)	Attention gating
Basal ganglia	Beta (13–30 Hz)	High-gamma	Motor control
Prefrontal	Delta (1–4 Hz)	Theta	Working memory

2.3 Role in SCPN¶

PAC measures hierarchical coupling between SCPN layers with different timescales. For example:

Layer 4 (Synchronisation, τ=2s) phase → Layer 2 (Neurochemical, τ=25ms) amplitude: slow oscillatory phase gates fast neurochemical activity
Layer 6 (Planetary, τ=24h) phase → Layer 5 (Psychoemotional, τ=1s) amplitude: circadian modulation of emotional reactivity

The pac_max metric (maximum PAC across all layer pairs) is consumed by supervisor/policy_rules.py for regime classification.

2.4 Comparison with Other Coupling Measures¶

Measure	What it captures	Frequency specificity	Directionality
PLV	Phase-phase coupling	Same frequency	Symmetric
PAC (MI)	Phase-amplitude coupling	Cross-frequency	Asymmetric
Power correlation	Amplitude-amplitude	Same or cross	Symmetric
Transfer entropy	Information flow	Any	Directed
Granger causality	Linear prediction	Any	Directed

PAC is unique in capturing cross-frequency interactions. PLV only works at the same frequency; transfer entropy is general but computationally expensive and harder to interpret.

2.5 Cross-Frequency Coupling Zoo¶

PAC is one member of a family of cross-frequency coupling (CFC) measures:

CFC type	Phase of	Amplitude of	Measure
Phase-amplitude (PAC)	Low freq	High freq	MI
Phase-phase (PPC)	Low freq	High freq	PLV at harmonics
Amplitude-amplitude (AAC)	Low freq	High freq	Correlation
Phase-frequency (PFC)	Low freq	—	Instantaneous freq vs phase

PAC is the most studied because it has the clearest functional interpretation: slow rhythms "gate" fast activity, enabling multiplexed information routing.

2.6 Methodological Considerations¶

Filtering artefacts. Band-pass filtering to extract low-frequency phase and high-frequency amplitude can introduce spurious PAC through filter edge effects. The Hilbert transform used in SPO's market.extract_phase() avoids some of these issues but is sensitive to broadband signals.

Volume conduction. In EEG/MEG, PAC measured between nearby channels may reflect volume conduction rather than genuine coupling. Source localisation (beamforming, LCMV) before PAC computation is recommended.

Harmonic confounds. A non-sinusoidal waveform (e.g., sharp-wave ripple) produces harmonics that can be mistaken for PAC. The "waveform shape" test (Cole & Voytek 2017) should be applied before interpreting PAC as genuine cross-frequency interaction.

2.7 PAC in Disease¶

Condition	PAC change	Brain region	Clinical relevance
Parkinson's disease	β→γ PAC elevated	Basal ganglia	Motor symptom correlate
Epilepsy	θ→γ PAC during seizure	Hippocampus	Seizure marker
Alzheimer's	θ→γ PAC reduced	Cortex	Memory impairment
Schizophrenia	θ→γ PAC disrupted	Prefrontal	Working memory deficit

These findings motivate PAC as a biomarker for closed-loop neurostimulation — the SCPN framework could modulate layer coupling based on real-time PAC.

3. Pipeline Position¶

┌──────────────┐     ┌──────────────────┐     ┌──────────────┐
│ Phase history │────→│ modulation_index │────→│ MI ∈ [0, 1] │
│ θ_low (T,)   │     │ (Tort et al.     │     └──────────────┘
│              │     │  2010)           │
│ Amplitude    │────→│                  │     ┌──────────────┐
│ A_high (T,)  │     └──────────────────┘     │ pac_gate     │
└──────────────┘                              │ MI > thresh? │
                     ┌──────────────────┐     └──────────────┘
┌──────────────┐     │ pac_matrix       │
│ θ (T, N)     │────→│ N×N MI matrix   │────→│ supervisor   │
│ A (T, N)     │────→│                  │     │ policy_rules │
└──────────────┘     └──────────────────┘     └──────────────┘

Inputs: - theta_low (T,) — low-frequency phase time series - amp_high (T,) — high-frequency amplitude time series - n_bins — number of phase bins (default 18)

Outputs: - modulation_index → float \(\in [0, 1]\) - pac_matrix → (N, N) MI matrix - pac_gate → bool

4. Features¶

4.1 Three Functions¶

Function	Purpose	Input	Output
`modulation_index`	Pairwise PAC	`(T,)` + `(T,)`	float MI
`pac_matrix`	All-pairs PAC	`(T,N)` + `(T,N)`	`(N,N)`
`pac_gate`	Binary decision	float MI	bool

4.2 Accelerator Output Contract¶

modulation_index delegates to spo_kernel.pac_modulation_index when available. The Rust implementation uses sorted-bin counting instead of np.digitize, avoiding branch mispredictions. Optional Rust, Go, Julia, and Mojo returns use the same output contract before public publication: modulation-index scalars must be finite real values in [0, 1], and flattened PAC-matrix payloads must contain exactly N*N finite real values in [0, 1]. Malformed optional backend output raises an error instead of being clipped into synthetic PAC evidence.

Public phase/amplitude vectors and histories reject boolean, complex, and numeric-string aliases before float conversion. Direct Go, Julia, and Mojo inputs use the same real-numeric contract. The Julia PAC-matrix bridge keeps raw return dtype intact until the shared N*N unit-interval validator runs, so stringified matrix entries cannot be laundered into PAC evidence. Mojo stdout remains an explicit strictly parsed text protocol.

4.3 Robustness¶

Empty arrays → MI = 0.0
n_bins < 2 or boolean bin aliases → ValueError
n_bins = 1 → rejected before the \(\log(1) = 0\) denominator
All-zero amplitude → MI = 0.0
Mismatched lengths → uses min(len(theta), len(amp))

5. Usage Examples¶

5.1 Basic Modulation Index¶

import numpy as np
from scpn_phase_orchestrator.upde.pac import modulation_index

T = 1000
t = np.arange(T) * 0.001  # 1ms steps

# Low-frequency phase (10 Hz)
theta_low = (2 * np.pi * 10 * t) % (2 * np.pi)

# High-frequency amplitude modulated by low-freq phase
amp_high = 1.0 + 0.5 * np.cos(theta_low)  # coupled!
amp_high += np.random.default_rng(42).normal(0, 0.1, T)  # noise

mi = modulation_index(theta_low, amp_high, n_bins=18)
print(f"MI = {mi:.4f}")  # expect MI > 0.3 (genuine coupling)

5.2 No Coupling (Null Case)¶

# Independent signals — MI should be near 0
rng = np.random.default_rng(0)
theta_random = rng.uniform(0, 2*np.pi, T)
amp_random = np.abs(rng.normal(0, 1, T))

mi_null = modulation_index(theta_random, amp_random)
print(f"Null MI = {mi_null:.4f}")  # expect < 0.05

5.3 PAC Matrix Across Layers¶

from scpn_phase_orchestrator.upde.pac import pac_matrix

N = 4  # 4 oscillators
T = 2000
phases = rng.uniform(0, 2*np.pi, (T, N))
amplitudes = np.abs(rng.normal(1, 0.3, (T, N)))

# Add coupling: phase[0] modulates amplitude[2]
amplitudes[:, 2] += 0.3 * np.cos(phases[:, 0])

pac = pac_matrix(phases, amplitudes, n_bins=18)
print("PAC matrix:")
print(np.round(pac, 3))
# Entry [0, 2] should be elevated

5.4 PAC Gate for Policy¶

from scpn_phase_orchestrator.upde.pac import pac_gate

mi = modulation_index(theta_low, amp_high)
if pac_gate(mi, threshold=0.3):
    print("Significant PAC detected — cross-frequency coupling active")
else:
    print("No significant PAC")

5.5 Bin Count Sensitivity¶

for B in [6, 9, 18, 36]:
    mi = modulation_index(theta_low, amp_high, n_bins=B)
    print(f"B={B:2d}: MI={mi:.4f}")
# MI should be consistent (± 0.05) across reasonable B values

5.6 Surrogate Testing¶

# Compute significance via phase shuffling
mi_real = modulation_index(theta_low, amp_high)
mi_surrogates = []
for i in range(200):
    # Shuffle amplitude to destroy PAC while preserving spectra
    shift = rng.integers(50, T - 50)
    amp_shuffled = np.roll(amp_high, shift)
    mi_surrogates.append(modulation_index(theta_low, amp_shuffled))

mi_surr = np.array(mi_surrogates)
z_score = (mi_real - mi_surr.mean()) / mi_surr.std()
p_value = np.mean(mi_surr >= mi_real)
print(f"MI = {mi_real:.4f}, z = {z_score:.2f}, p = {p_value:.4f}")

5.7 SCPN Cross-Layer PAC¶

# Layer 4 (sync) phase → Layer 2 (neurochemical) amplitude
from scpn_phase_orchestrator.upde.engine import UPDEEngine
from scpn_phase_orchestrator.upde.order_params import compute_order_parameter

# (Assume simulation has been run, producing phase trajectories per layer)
# layer4_phases = phase_trajectory[:, layer4_indices]  # (T, N_layer4)
# layer2_amplitudes = amplitude_trajectory[:, layer2_indices]

# Cross-layer PAC: does Layer 4 phase gate Layer 2 amplitude?
# pac_cross = pac_matrix(layer4_phases, layer2_amplitudes)

5.8 Controlled Coupling Strength¶

# Generate signals with known MI for validation
coupling_strengths = [0.0, 0.1, 0.3, 0.5, 0.8, 1.0]
for c in coupling_strengths:
    theta = (2 * np.pi * 10 * t) % (2 * np.pi)
    amp = 1.0 + c * np.cos(theta) + rng.normal(0, 0.1, T)
    mi = modulation_index(theta, np.abs(amp))
    print(f"Coupling={c:.1f}: MI={mi:.4f}")
# MI should increase monotonically with coupling strength

5.9 PAC Matrix Visualisation¶

import matplotlib.pyplot as plt

pac = pac_matrix(phases_history, amplitudes_history)
fig, ax = plt.subplots()
im = ax.imshow(pac, cmap="hot", vmin=0, vmax=1)
ax.set_xlabel("Amplitude oscillator j")
ax.set_ylabel("Phase oscillator i")
ax.set_title("Phase-Amplitude Coupling Matrix")
plt.colorbar(im, label="MI")
plt.savefig("pac_matrix.png", dpi=150)

5.10 Time-Resolved PAC¶

# Sliding-window PAC to track coupling over time
window = 500
mi_timeseries = []
for start in range(0, T - window, 50):
    end = start + window
    mi = modulation_index(theta_low[start:end], amp_high[start:end])
    mi_timeseries.append(mi)
# mi_timeseries shows how PAC evolves — useful for detecting
# coupling onset/offset in SCPN regime transitions

6. Technical Reference¶

6.1 Module API¶

pac ¶

Tort 2010 phase-amplitude coupling with 5-backend fallback chain.

Follows feedback_module_standard_attnres.md:

modulation_index — scalar Tort 2010 MI on a single (θ_low, a_high) pair of time series.
pac_matrix — (N, N) pairwise MI matrix over N oscillator phase / amplitude channels.
pac_gate — pure-Python boolean gate on an MI value (no backend dispatch needed — trivial comparison).

All compute kernels are available in Rust, Mojo, Julia, Go, Python. AVAILABLE_BACKENDS reports detected backends in canonical fallback order, while ACTIVE_BACKEND is selected by a small import-time hot-path probe so slow external wrappers do not displace the faster local path.

Functions:¶

modulation_index ¶

modulation_index(
    theta_low: FloatArray,
    amp_high: FloatArray,
    n_bins: int = 18,
) -> float

Phase-amplitude coupling via Tort et al. 2010, J. Neurophysiol.

Bins amplitude by phase, computes KL divergence from uniform, returns the modulation index normalised to [0, 1] by log(n_bins).

Parameters¶

theta_low : FloatArray Low-frequency driver phase in radians, shape (T,). amp_high : FloatArray High-frequency amplitude envelope, shape (T,). n_bins : int Number of phase bins used for the modulation-index histogram.

Returns¶

float The Tort modulation index of phase-amplitude coupling.

Raises¶

ValueError If n_bins is not a positive integer or inputs mismatch.

Source code in src/scpn_phase_orchestrator/upde/pac.py

def modulation_index(
    theta_low: FloatArray, amp_high: FloatArray, n_bins: int = 18
) -> float:
    """Phase-amplitude coupling via Tort et al. 2010, J. Neurophysiol.

    Bins amplitude by phase, computes KL divergence from uniform,
    returns the modulation index normalised to ``[0, 1]`` by
    ``log(n_bins)``.

    Parameters
    ----------
    theta_low : FloatArray
        Low-frequency driver phase in radians, shape ``(T,)``.
    amp_high : FloatArray
        High-frequency amplitude envelope, shape ``(T,)``.
    n_bins : int
        Number of phase bins used for the modulation-index histogram.

    Returns
    -------
    float
        The Tort modulation index of phase-amplitude coupling.

    Raises
    ------
    ValueError
        If ``n_bins`` is not a positive integer or inputs mismatch.
    """
    n_bins = _validate_n_bins(n_bins)
    theta_low = _validate_signal("theta_low", theta_low)
    amp_high = _validate_signal("amp_high", amp_high)
    if np.any(amp_high < 0.0):
        raise ValueError("amp_high must contain non-negative amplitudes")
    if theta_low.size == 0 or amp_high.size == 0:
        return 0.0

    backend_fn = _dispatch("modulation_index")
    if backend_fn is not None:
        fn = cast(
            "Callable[[FloatArray, FloatArray, int], float]",
            backend_fn,
        )
        return _validate_mi_value(
            "modulation_index backend",
            fn(
                np.ascontiguousarray(theta_low, dtype=np.float64),
                np.ascontiguousarray(amp_high, dtype=np.float64),
                n_bins,
            ),
        )

    return _validate_mi_value(
        "modulation_index",
        _modulation_index_python(theta_low, amp_high, n_bins),
    )

pac_matrix ¶

pac_matrix(
    phases_history: FloatArray,
    amplitudes_history: FloatArray,
    n_bins: int = 18,
) -> FloatArray

Return the (N, N) PAC matrix [i, j] = MI(phase_i, amplitude_j).

Parameters¶

phases_history : FloatArray (T, N) phase time series. amplitudes_history : FloatArray (T, N) amplitude time series. n_bins : int number of phase bins.

Returns¶

FloatArray FloatArray The (N, N) phase-amplitude coupling matrix.

Raises¶

ValueError If n_bins is not positive or the histories have mismatched shapes.

Source code in src/scpn_phase_orchestrator/upde/pac.py

def pac_matrix(
    phases_history: FloatArray,
    amplitudes_history: FloatArray,
    n_bins: int = 18,
) -> FloatArray:
    """Return the ``(N, N)`` PAC matrix ``[i, j] = MI(phase_i, amplitude_j)``.

    Parameters
    ----------
    phases_history : FloatArray
        ``(T, N)`` phase time series.
    amplitudes_history : FloatArray
        ``(T, N)`` amplitude time series.
    n_bins : int
        number of phase bins.

    Returns
    -------
    FloatArray
        FloatArray The ``(N, N)`` phase-amplitude coupling matrix.

    Raises
    ------
    ValueError
        If ``n_bins`` is not positive or the histories have mismatched shapes.
    """
    n_bins = _validate_n_bins(n_bins)
    phases_history = _validate_history("phases_history", phases_history)
    amplitudes_history = _validate_history("amplitudes_history", amplitudes_history)
    if np.any(amplitudes_history < 0.0):
        raise ValueError("amplitudes_history must contain non-negative amplitudes")
    t, n = phases_history.shape
    if t == 0 or n == 0:
        return np.zeros((n, n), dtype=np.float64)
    if amplitudes_history.shape != (t, n):
        raise ValueError("phases and amplitudes must have the same shape")

    backend_fn = _dispatch("pac_matrix")
    if backend_fn is not None:
        fn = cast(
            ("Callable[[FloatArray, FloatArray, int, int, int], FloatArray]"),
            backend_fn,
        )
        flat = fn(
            np.ascontiguousarray(phases_history.ravel(order="C"), dtype=np.float64),
            np.ascontiguousarray(amplitudes_history.ravel(order="C"), dtype=np.float64),
            t,
            n,
            n_bins,
        )
        matrix = validate_pac_matrix_output(flat, n=n)
        return matrix.reshape((n, n), order="C")

    result = np.zeros((n, n), dtype=np.float64)
    for i in range(n):
        for j in range(n):
            result[i, j] = modulation_index(
                phases_history[:, i], amplitudes_history[:, j], n_bins
            )
    return result

pac_gate ¶

pac_gate(pac_value: float, threshold: float = 0.3) -> bool

Binary gate: True when PAC exceeds threshold.

Pure-Python helper; no dispatcher — the comparison is trivial.

Parameters¶

pac_value : float A phase-amplitude coupling value. threshold : float Decision threshold.

Returns¶

bool True when the PAC value exceeds the threshold.

Source code in src/scpn_phase_orchestrator/upde/pac.py

def pac_gate(pac_value: float, threshold: float = 0.3) -> bool:
    """Binary gate: ``True`` when PAC exceeds ``threshold``.

    Pure-Python helper; no dispatcher — the comparison is trivial.

    Parameters
    ----------
    pac_value : float
        A phase-amplitude coupling value.
    threshold : float
        Decision threshold.

    Returns
    -------
    bool
        ``True`` when the PAC value exceeds the threshold.
    """
    pac_value = _validate_finite_real("pac_value", pac_value)
    threshold = _validate_finite_real("threshold", threshold)
    return pac_value >= threshold

6.2 Function Signatures¶

modulation_index(theta_low: NDArray, amp_high: NDArray, n_bins: int = 18) → float

Parameter	Type	Default	Description
`theta_low`	`NDArray` (T,)	—	Low-frequency phase
`amp_high`	`NDArray` (T,)	—	High-frequency amplitude
`n_bins`	`int`	`18`	Phase bins
Returns	`float`		MI \(\in [0, 1]\)

pac_matrix(phases_history: NDArray, amplitudes_history: NDArray, n_bins: int = 18) → NDArray

Parameter	Type	Description
`phases_history`	`(T, N)`	Phase time series
`amplitudes_history`	`(T, N)`	Amplitude time series
Returns	`(N, N)`	MI matrix

pac_gate(pac_value: float, threshold: float = 0.3) → bool

7. Performance Benchmarks¶

7.1 modulation_index¶

T	n_bins	Python (µs)	Rust (µs)	Speedup
500	18	85	12	7.1x
5000	18	420	55	7.6x
50000	18	4200	480	8.8x

7.2 pac_matrix¶

\(O(N^2)\) calls to modulation_index:

N	T	Python (ms)	Rust (ms)	Speedup
4	1000	1.4	0.2	7x
16	1000	22	3.1	7.1x
64	1000	350	48	7.3x

7.3 Complexity¶

Function	Time	Space
`modulation_index`	\(O(T + B)\)	\(O(T + B)\)
`pac_matrix`	\(O(N^2 \cdot (T + B))\)	\(O(N^2 + T + B)\)
`pac_gate`	\(O(1)\)	\(O(1)\)

7.4 Surrogate Testing Cost¶

200 surrogates × MI computation:

T	Python (ms)	Rust (ms)
1000	17	2.4
5000	84	11
50000	840	96

7.5 pac_matrix Scaling¶

The matrix computation is \(O(N^2)\) MI evaluations. For large \(N\), this becomes the bottleneck:

N	T=1000, Python	T=1000, Rust
4	1.4 ms	0.2 ms
8	5.4 ms	0.8 ms
16	22 ms	3.1 ms
32	87 ms	12 ms
64	350 ms	48 ms

For \(N > 32\), consider computing PAC only for specific layer pairs rather than the full matrix.

7.6 Recommended Settings¶

Use case	T	n_bins	Surrogates	Expected time (Rust)
Quick check	500	18	0	0.01 ms
Standard analysis	2000	18	200	4 ms
Publication	10000	18	500	25 ms
Full matrix (N=16)	2000	18	200 per pair	~1s

7.7 Memory¶

MI computation is \(O(T + B)\) — the phase binning creates \(B\) accumulators and a bin-index array of length \(T\). For \(T = 50000\) and \(B = 18\): ~400 KB. Negligible compared to input data.

8. Citations¶

Tort A.B.L., Komorowski R., Eichenbaum H., Kopell N. (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. Journal of Neurophysiology 104(2):1195–1210. doi:10.1152/jn.00106.2010
Canolty R.T., Edwards E., Dalal S.S., et al. (2006). High gamma power is phase-locked to theta oscillations in human neocortex. Science 313(5793):1626–1628. doi:10.1126/science.1128115
Jensen O., Colgin L.L. (2007). Cross-frequency coupling between neuronal oscillations. Trends in Cognitive Sciences 11(7):267–269. doi:10.1016/j.tics.2007.05.003
Aru J., Aru J., Priesemann V., et al. (2015). Untangling cross- frequency coupling in neuroscience. Current Opinion in Neurobiology 31:51–61. doi:10.1016/j.conb.2014.08.002
Hülsemann M.J., Naumann E., Rasch B. (2019). Quantification of phase-amplitude coupling in neuronal oscillations: comparison of phase-locking value, mean vector length, modulation index, and generalized-linear-modeling-cross-frequency-coupling. Frontiers in Neuroscience 13:573. doi:10.3389/fnins.2019.00573
Cole S.R., Voytek B. (2017). Brain oscillations and the importance of waveform shape. Trends in Cognitive Sciences 21(2):137–149. doi:10.1016/j.tics.2016.12.008
Hyafil A., Giraud A.-L., Fontolan L., Gutkin B. (2015). Neural cross-frequency coupling: connecting architectures, mechanisms, and functions. Trends in Neurosciences 38(11):725–740. doi:10.1016/j.tins.2015.09.001

Test Coverage¶

tests/test_pac.py — 20 tests: MI bounds, non-negative amplitude enforcement, vector/history validation, dispatcher fallback, backend payload validation, pac_matrix shape/range checks, and pac_gate threshold semantics.
tests/test_pac_parity.py — 20 tests: Python fallback contracts, malformed inputs, backend fallback payload checks, Rust parity when available, and matrix layout guards.
tests/test_pac_backends.py — 112 tests: direct Rust, Julia, Go, and Mojo boundary contracts, runtime-loading guards, parity checks, Mojo text-protocol cardinality, and cross-backend consistency.
tests/test_pac_stability.py — 7 slow stability/property tests for bounded MI, modulation-depth ordering, uncorrelated noise floor, and matrix invariants.
tests/test_pac_public_backend_output_contracts.py — 6 public dispatcher tests covering optional backend output rejection, valid output publication, loader-failure fallback, and timing-probe unavailability.

Total collected PAC tests: 165.

Multi-backend fallback chain¶

Since the 2026-04-17 migration to the AttnRes-level module standard, pac ships with five language-backed implementations for every kernel (modulation_index and pac_matrix). The dispatcher resolves the fastest available backend at import time and exposes the choice as ACTIVE_BACKEND / AVAILABLE_BACKENDS.

Position	Backend	Build
1	Rust	`maturin develop -m spo-kernel/crates/spo-ffi/Cargo.toml --release`
2	Mojo	`mojo build mojo/pac.mojo -o mojo/pac_mojo -Xlinker -lm`
3	Julia	`juliacall` + `julia/pac.jl`
4	Go	`cd go && go build -buildmode=c-shared -o libpac.so pac.go`
5	Python	always present

Parity (measured vs NumPy reference)¶

Backend	`modulation_index`
Rust	6.59e-17 (bit-exact)
Mojo	3.70e-12 (text-protocol budget)
Julia	6.94e-18 (bit-exact)
Go	6.94e-17 (bit-exact)
Python	0 (reference)

Tests in tests/test_pac_backends.py enforce atol = 1e-12 for Rust/Julia/Go and atol = 1e-10 for Mojo (the log-sum-over-bins step amplifies the 17-digit text-protocol floor).

Measured benchmark¶

Output from PYTHONPATH=src python benchmarks/pac_benchmark.py --sizes 200 1000 5000 --calls 100:

N	Rust	Mojo	Julia	Go	Python
200	0.078 ms	111.378 ms	80.976 ms	0.241 ms	0.332 ms
1000	0.010 ms	99.913 ms	0.032 ms	0.028 ms	0.328 ms
5000	0.035 ms	111.358 ms	0.053 ms	0.045 ms	0.773 ms

Rust dominates; Julia/Go converge to within 1.3× of Rust after JIT warm-up at N ≥ 1000. NumPy fallback stays within 22× of Rust — usable for one-off analysis but not for a monitor loop running on every simulation step. Mojo's subprocess floor rules it out of the hot loop until Mojo 0.27+.

Source¶

Python dispatcher: src/scpn_phase_orchestrator/upde/pac.py
Python bridges: upde/_pac_julia.py, upde/_pac_go.py, upde/_pac_mojo.py
Rust: spo-kernel/crates/spo-engine/src/pac.rs
Julia: julia/pac.jl
Go: go/pac.go
Mojo: mojo/pac.mojo
FFI: spo-kernel/crates/spo-ffi/src/lib.rs (pac_modulation_index, pac_matrix_compute)
Tests: tests/test_pac.py, tests/test_pac_parity.py, tests/test_pac_backends.py, tests/test_pac_stability.py
Benchmark: benchmarks/pac_benchmark.py