Biphasic Hierarchical Reorganization and Riemannian Connectomic Plasticity: A Longitudinal Precision Functional Mapping Study of Psilocybin

1. Introduction: From Euclidean Artifacts to Riemannian Manifolds

Contemporary psychiatry is experiencing a renaissance in serotonergic psychedelic research, driven by the imperative to treat refractory affective disorders. However, the computational frameworks underpinning this movement remain largely constrained by Euclidean statistical methodologies that assume a flat, well-behaved error variance. Such assumptions undergo geometric collapse under the massive entropic impact of $5\text{-HT}_{2\text{A}}$ receptor activation, leading to methodological myopia (Arsigny et al., 2006) and frequent topological false negatives (Nichols et al., 2017). Traditional linear tools are algebraically equivalent to using a straight ruler to measure distance on the surface of a sphere, an approach that artificially inflates connectome variance and masks subacute traits of plasticity, wherein the true magnitude of structural brain alterations is frequently obscured by computational noise.

The central objective of this work is to diverge from this tradition by treating the human connectome as a rigid mathematical topology resident on a Riemannian manifold. Utilizing a dense longitudinal precision functional mapping (PFM) dataset of $N = 101$ fMRI sessions, I aimed to precisely map the kinematics of psilocybin-induced neuroplasticity across three strict temporal windows: prior homeostasis (Baseline), the entropic storm (Acute Phase), and the structural settling (Persistent Phase).

To ensure inferential transparency and anatomical fidelity, I employed Glasser/CABNP functional parcellation. The computational pipeline implements a mandatory medial wall structural purge executed at the high-density grayordinate scale, excluding 76 non-tissue vertices from the raw surface input substrate (reducing the local coordinate grid from 91,282 to 91,206 grayordinates) to completely prevent midline signal pollution. Downstream spatial projection onto the subcortical-inclusive CABNP functional configuration maps these filtered vectors directly onto its optimized, invariant 717-node partition. This routine enables the precise tracking of cortical hubs along their exact stereotaxic axes while maintaining total dimensional parity across all longitudinal testing lines. By enforcing Symmetric Positive-Definite (SPD) manifold integrity via Ledoit-Wolf shrinkage, this study unearths the pure biological signal of long-term cortical remodeling that survives weeks after the complete metabolization of the compound.

To prevent the temporal mapping of this neuroplastic window from collapsing under the weight of data instability, it was imperative to rectify a foundational flaw in classical neuroimaging. Functional matrices under the influence of psychedelics inject a surge of heteroscedastic variance into the system. Therefore, I employed Riemannian geometry coupled with the Ledoit-Wolf shrinkage estimator and robust logarithmic projections (HC3). This mathematical shielding isolated the pure biological signal, revealing structural dynamics obscured by previous methods.

The Four Pillars of Translational Discovery

By aligning computational rigor with clinical outcomes, this study establishes a translational overview grounded in four fundamental pillars:

2. Methods: Precision Functional Mapping Architecture

To evaluate potential psilocybin-induced topological modifications, the analytical framework must address a well-documented challenge in neuroimaging: the high-dimensionality problem. When estimating functional brain networks using extensive spatial parcellations over discrete time intervals, the resulting empirical covariance matrices can become ill-conditioned (Varoquaux & Craddock, 2013). Spatial parcellations that omit subcortical structures or fail to adequately segregate inter-hemispheric boundaries may obscure genuine topological properties (Glasser et al., 2016; Ji et al., 2019).

If the empirical covariance matrix lacks full-rank status or contains non-positive eigenvalues, non-Euclidean geometric frameworks are susceptible to the tensor-swelling effect. Ensuring symmetric positive-definite (SPD) manifold integrity offers a mathematically consistent basis for tracking long-term geodesic distances across temporal windows (Arsigny et al., 2006). To address these limitations, I implemented a structured pipeline incorporating anatomical spatial anchoring, covariance regularization, and robust variance estimation.

Cohort Matrix & Longitudinal Architecture

Data were obtained from a dense-sampling precision functional mapping framework designed to track the dynamic transitions associated with psilocybin administration:

• Cohort Matrix: The sample consisted of 7 healthy adult volunteers ($M = 34.1\text{ years}$, $SD = 9.8$; $n = 3\text{ females}$; 6 participants identified as Caucasian).
• Longitudinal Architecture: The imaging protocol spanning 11 separate longitudinal subject lines captured three distinct temporal windows: baseline homeostasis (Baseline), the acute drug effect window (60–90 min post-ingestion; Acute Phase), and a subacute settling window (up to 2 weeks post-exposure; Persistent Phase).
• Replication Protocol: To evaluate the consistency of these network trajectories, a subset of 4 individuals (P1, P3, P4, P5) underwent a replication protocol following a washout period of greater than 6 months. In alignment with precision functional mapping conventions, these replication sequences were modeled as distinct longitudinal streams, yielding a total dataset of 101 dense imaging sessions.

Multi-Echo EPI Acquisition & NORDIC Denoising

Imaging data were acquired on a 3T Siemens Prisma system:

• Sequence: Multi-echo echo-planar imaging (EPI; five TEs: 14.20, 38.93, 63.66, 88.39, 113.12 ms) acquired with a repetition time (TR) of 1761 ms, flip angle = 68°, and in-plane acceleration (iPAT/GRAPPA) = 2, optimized for absolute physical signal separation.
• Normalization: Blood-oxygen-level-dependent (BOLD) time-series were subjected to thermal and non-biological noise suppression via the NORDIC denoising algorithm, followed by temporal Z-score standardization (zero mean and unit variance) for amplitude calibration. For runs truncated due to acute session termination (e.g., P6 ses-5 REST3; P4R ses-2 REST3), the 3 missing trailing noise calibration frames were appended from adjacent, identical resting-state runs within the same subject matrix. This adjustment was strictly restricted to non-signal calibration volumes; the 513-frame internal temporal autocorrelation structure and phase-amplitude relationships of the primary BOLD time-series were left entirely unaltered, shielding the underlying Riemannian manifold geometry from interpolation artifacts (Nichols et al., 2017).

Parcellation, Medial Wall Purge, and Spatial Anchoring

To minimize signal cross-contamination, time-series were extracted by mapping a high-density 91,206 grayordinate substrate onto a validated global cortico-subcortical parcellation consisting of an optimized 717-node partition derived from the CABNP atlas (Ji et al., 2019). The functional signal was extracted by applying a cortical ribbon-constrained interpolation method.

To limit partial volume effects along the brain midline, a structural mask was applied to exclude exactly 76 non-tissue vertices corresponding to the medial wall boundary (Glasser et al., 2016). This spatial filtering ensures that the resulting nodal coordinates maintain consistent stereotaxic references across sessions.

Riemannian Metric & Covariance Regularization

Direct empirical estimation of a covariance network containing 717 nodes ($p = 717$) from finite temporal samples ($n$) typically yields ill-conditioned or singular matrices (Varoquaux & Craddock, 2013). Measuring distances between brain states using standard linear metrics can introduce numerical instability and inflate connectome variance, an effect known as tensor-swelling (Arsigny et al., 2006).

To stabilize the network matrices against high-frequency fluctuations, I implemented the Ledoit-Wolf shrinkage estimator (Ledoit & Wolf, 2004). Across all conditions, this regularization maintained full-rank integrity, satisfying the non-singular requirement for Riemannian operations with an absolute session-wise eigenvalue floor of $\lambda_{\min} > 0$. Group consensus condition matrices converged at a higher boundary threshold ($\lambda_{\min} = 1.042 \times 10^{-4}$). This preconditioning provides a stable foundation for computing geodesic distances using the Affine-Invariant Riemannian Metric (AIRM) across experimental windows (Pennec et al., 2006).

3. Results: Topological Transformations Across Temporal Windows

Riemannian Integrity & Discovery of Biphasic Hierarchical Reorganization

The regularization architecture successfully secured full-rank status and symmetric positive-definite (SPD) manifold integrity across all 101 sessions ($\lambda_{\min} > 0$). Group-level consensus condition matrices demonstrated stable convergence boundaries: baseline surface template recorded a minimum eigenvalue of $\lambda_{\min} = 1.042 \times 10^{-4}$, the acute contractive matrix yielded $\lambda_{\min} = 1.037 \times 10^{-4}$, and the subacute persistent phase settled at $\lambda_{\min} = 1.042 \times 10^{-4}$.

Analysis of low-dimensional spectral embeddings tracked variations along the principal functional gradient across time. The distance between transmodal and unimodal poles, which stood at a mean of 0.002707 at Baseline, contracted to 0.000841 during the Acute Phase ($\Delta = -68.95\%$), reflecting a transient compression of the functional hierarchy. During the post-acute Persistent Phase, the network layout exhibited a subacute expansion, reaching a polar distance mean of 0.002947 ($\Delta = +250.66\%$ relative to the acute state).

Global Disintegration and the Acute Connective Avalanche

During peak acute exposure, the local segregation index of the Default Mode Network (DMN) underwent a precipitous drop from a Baseline mean of .4608 to an Acute mean of .2547. Concurrently, the global mean participation coefficient increased from .1910 to .2816. Global network deformation (RiemannianPlasticity) surged from a baseline mean of 34.0016 to an acute peak of 53.0813.

In the persistent follow-up window, DMN segregation recovered to a mean of 0.4800, while global deformation settled at 43.7956, characterized by a drastic tightening of individual session variance ($SD = 4.0266$, compared to baseline $SD = 10.5015$).

The Phenomenological Interaction Space (MEQ-30)

I tested the relationship between subjective mystical scores and objective connectomic remodeling. The cohort exhibited broad subjective variance on the MEQ-30 ($M = 3.0697$, $SD = 1.2623$, range 1.1333 – 4.4000).

Longitudinal GEE modeling revealed that the long-term main effect of the mean-centered MEQ-30 score was statistically non-significant ($b = 0.0315$, $p = .513$). However, a highly significant trait-by-state interaction emerged during the peak acute window: acute GradientShift displacement was strongly predicted by the interaction between temporal state and centered psychometrics ($b = 0.1550$, $p = .008$). In the persistent window, this predictive value dissipated ($b = 0.0341$, $p = .137$).

Alternative Regression Frameworks & Heteroscedastic Robustness

To account for non-constant variance without data-trimming, I compared standard OLS (Model D.2) against a regularized Gamma Generalized Linear Model with log-link and Huber-White HC3 sandwich covariance estimation (Model D.1).

Both models converged on identical tracking behavior: acute perturbation was highly significant ($p < .001$), while the persistent phase estimator crossed the zero-effect boundary ($p = .963$ in OLS, $p = .981$ in Gamma GLM).

Subnetwork Trajectories & Transmodal Asymmetry (Depression vs. Anxiety)

Evaluating 12 transmodal control nodes (DMN/FPN mood regulation clusters) versus 12 salience-limbic nodes (insula/amygdala hypervigilance complexes) revealed profound structural divergence:

• Transmodal cortical hub Default-66_L-Ctx exhibited a localized displacement vector of 0.2133 ($\Delta = 0.2963$ Baseline → Acute).
• Salience hub Cingulo-Opercular-42_L-Ctx reached only 0.1255 ($\Delta = 0.1901$ Acute → Persistent).

A non-parametric Mann-Whitney U test decisively rejected isotropic spatial distribution: $U = 139.00$, $p < .001$, rank-biserial correlation $r = -.9306$. Transmodal depression-related circuits absorb significantly greater coordinate malleability than salience or limbic anxiety pathways.

4. Discussion: Translating Connectomic Geometry into Translational Biology

Having established the topological and statistical basis of psilocybin-induced neuronal displacement, the imperative now is to translate tensors and Riemannian manifolds into translational biology. Decoding the acute compression and subsequent persistent hyperexpansion of cortical networks demands coupling this architectural deformation to the primary pharmacodynamics of $5\text{-HT}_{2\text{A}}$ receptors and defining its impact on patient prognosis.

Neuronal Displacement and REBUS Validation

The coordinate displacement vectors identified along the principal tracking axis point to clear adjustments within macroscale cortical hierarchies. The observed contraction in polar distance (0.002707 → 0.000841, $\Delta = -68.95\%$) indicates a transient reduction in macroscopic hierarchical differentiation rather than an unconstrained artifact.

Under resting-state conditions, the cerebral cortex maintains organized topological boundaries separating primary sensory streams from transmodal higher-order networks (Margulies et al., 2016). In affective disorders, these higher-order systems become overly rigid, locking patients into depressive rumination. The acute tracking data demonstrates a temporary mitigation of these structural boundaries: DMN segregation drops from .4608 to .2547, while mean participation rises from .1910 to .2816. This provides empirical verification of the REBUS (RElaxed Beliefs Under pSychedelics) framework (Carhart-Harris & Friston, 2019): $5\text{-HT}_{2\text{A}}$ agonism temporarily suppresses top-down priors, rendering underlying pathways exceptionally receptive to cognitive restructuring.

The Mystical Dogma: Phenomenological vs. Topological Reset

A prominent hypothesis in psychedelic medicine posits that clinical durability is causally driven by the acute subjective intensity of mystical experiences (Carhart-Harris et al., 2017). The longitudinal data from this study refines this view:

1. Long-term post-acute connectome remodeling is statistically orthogonal to the mean-centered MEQ-30 score ($p = .513$), supporting a standalone physiological reset mechanism.
2. In contrast, acute gradient displacement is strongly predicted by the trait-by-state interaction ($p = .008$).

Acute phenomenology directly reflects the immediate magnitude of topological destabilization, but long-term structural retention is governed by biological plasticity and individualized environmental interaction (Letheby, 2021).

Regional Susceptibility & The Subnetwork Asymmetry

The non-parametric discovery of transmodal asymmetry ($U = 139.00$, $p < .001$, $r = -.9306$) shatters the assumption of an isotropic, whole-brain reset. High-order transmodal assemblies (DMN/FPN; e.g., Default-66_L-Ctx, $M = 0.2133$) exhibit significantly greater structural malleability than lower-level salience or limbic structures (Cingulo-Opercular-42_L-Ctx, $M = 0.1255$).

This differential susceptibility aligns with cortical $5\text{-HT}_{2\text{A}}$ receptor density gradients, explaining why psilocybin produces rapid breakthroughs in ruminative depressive disorders while requiring specialized somatic protocols for panic or generalized anxiety.

5. Conclusion & Future Frontiers

In summary, regularized tracking across this 101-session precision functional mapping dataset establishes three fundamental kinetic phases:

1. Baseline Phase: Functional networks operate within typical modular boundaries, maintaining clear topological distinctions between primary sensory streams and transmodal control networks.
2. Acute Phase: Marked hierarchical gradient compression ($\Delta = -68.95\%$) and DMN desegregation (.4608 → .2547), creating high-entropy cross-system communication and relaxing top-down priors.
3. Persistent Phase: Subacute gradient expansion ($\Delta = +250.66\%$) and canalization into a tightened attractor state ($CV = 9.19\%$), retaining individualized malleability.

Methodological Limitations & Future Directions

• Sample Size vs. Precision: While 101 dense sessions provide unprecedented intra-subject statistical power, the sample is anchored to $N = 7$ unique participants. Validation in larger clinical cohorts is warranted.
• Multimodal hdEEG-fMRI Integration: Combining dense fMRI with high-density EEG will enable millisecond-resolved tracking of gradient collapse.
• Cross-Compound Comparison: Contrasting psilocybin against ketamine, LSD, and MDMA will determine whether transmodal asymmetry is unique to $5\text{-HT}_{2\text{A}}$ agonism or represents a general property of rapid psychoplastogens.

💻 Data and Code Availability

• Raw Neuroimaging Dataset: OpenNeuro Accession ID ds006072 (Psilocybin Precision Functional Mapping database).
• Computational Pipeline Source Code: The complete, version-controlled algorithmic pipeline (CABNP surface projection, Ledoit-Wolf shrinkage, Affine-Invariant Riemannian Metric tensors) is publicly available on GitLab:
  https://gitlab.com/lunalab/macroscale-cortical-gradient-disintegration-and-riemannian-connectomics-pipeline/-/tree/v.1.0.0?ref_type=tags
• Persistent Zenodo Archive: DOI: 10.5281/zenodo.22144006

📚 References