# paper.pdf ## Metadata - PDFFormatVersion=1.6 - Language=en-US - IsLinearized=true - IsAcroFormPresent=false - IsXFAPresent=false - IsCollectionPresent=false - IsSignaturesPresent=false - CreationDate=D:20260602190150+05'30' - Creator=Adobe InDesign 17.4 (Windows) - ModDate=D:20260602194822+05'30' - Producer=Adobe PDF Library 16.0.7 - Trapped=False - xmpmm:instanceid=uuid:3dcca18e-664e-4dda-b986-f6187cb1f53f - xmpmm:originaldocumentid=xmp.did:A33916180F2068118083ACD1B377E157 - xmpmm:documentid=xmp.id:477f070e-c9dd-db4a-8760-0da12a6391d0 - xmpmm:renditionclass=proof:pdf - xmpmm:derivedfrom=xmp.iid:3e67d309-bc91-594d-a07e-69d009957b81xmp.did:3e67d309-bc91-594d-a07e-69d009957b81xmp.did:A33916180F2068118083ACD1B377E157default - xmpmm:history=convertedfrom application/x-indesign to application/pdfAdobe InDesign 17.4 (Windows)/2026-06-02T19:01:50+05:30 - xmp:createdate=2026-06-02T19:01:50+05:30 - xmp:modifydate=2026-06-02T19:48:22+05:30 - xmp:metadatadate=2026-06-02T19:48:22+05:30 - xmp:creatortool=Adobe InDesign 17.4 (Windows) - xmprights:marked=False - photoshop:modifydate=2019-01-31T16:16:25 - photoshop:createdate=2019-01-24T11:43:46 - dc:format=application/pdf - pdf:producer=Adobe PDF Library 16.0.7 - pdf:trapped=False - extensisfontsense:slug=3023672065GraphikOpenType - PS3023672065--0Graphik-Bold1.0011904271029GraphikOpenType - PS1904271029--0Graphik-Regular1.0012428528031News Gothic MT OTOpenType - PS2428528031Monotype Typography0NewsGothicMTOT-Regular1.000543760319News Gothic MT OTOpenType - PS543760319Monotype Typography0NewsGothicMTOT-Bold1.0002785939051KnockoutOpenType - PS2785939051--0Knockout-68FullFeatherwt1.200423856041Graphik X CondOpenType - PS423856041--0GraphikXCond-Super1.0001526035435Graphik X CondOpenType - PS1526035435--0GraphikXCond-Bold1.0003093728306Graphik CondOpenType - PS3093728306--0GraphikCond-Black1.0011341118539Graphik CondOpenType - PS1341118539--0GraphikCond-Regular1.0011385893531GraphikOpenType - PS1385893531--0Graphik-Semibold1.0012986479693Graphik CompactOpenType - PS2986479693--0GraphikCompact-Regular1.0024223138361Graphik CompactOpenType - PS4223138361--0GraphikCompact-Semibold1.0023673047546Harding SevenOpenType - PS3673047546--0HardingSeven-Extrabold0.000103712566Harding SevenOpenType - PS103712566--0HardingSeven-Roman0.0002232452979Harding SevenOpenType - PS2232452979--0HardingSeven-Medium0.000631916388Harding SevenOpenType - PS631916388--0HardingSeven-Bold0.000157629484Publico BannerOpenType - PS157629484--0PublicoBanner-Black1.005 ## Contents ### Page 1 nature neurosciencehttps://doi.org/10.1038/s41593-026-02311-2ArticleHuman learning of noninvasive brain– computer interfaces via manifold geometryIn the format provided by the authors and uneditedSupplementary information ### Page 2 Human learning of noninvasive brain-computer interfaces via manifold geometrySupplementary informationThis file contains: Supplementary Methods Supplementary Table 1 Supplementary Figures 1 – 16 1 ### Page 3 Supplementary MethodsStimulus design and presentationWe used a custom video game task designed and programmed with C# for Unity (Version 2019.4.12f1). The video game was presented in the scanner on a screen at the back of the bore using a rear projector (1, 9201, 068-pixel resolution; 60-Hz refresh rate). Participants viewed the display through a mirror mounted on the MRI headcoil. The game was run through Unity on a Linux computer. fMRI scanner triggers were collected via PsychoPy (Version 2022.2.4) running in parallel with Unity on the Linux machine. Custom Python software connected PsychoPy and Unity via a TCP port to synchronize game events with scanner triggers. A message was sent from PsychoPy to Unity to start the game when the first scanner trigger was received, and messages were sent from Unity to PsychoPy at the start and end of each trial. In each trial of the game, participants were tasked with steering an avatar through an outdoor environment to a goal location (a flag). Landscape objects (shrubs, trees, and rocks) were randomly placed throughout the environment as landmarks and for visual engagement. Rocks and shrubs were lower than avatar knee-height, whereas trees were taller than avatar height. They were always rendered outside a wide buffer from the direct path between the initial and goal locations so they did not obstruct the path, but if participants strayed significantly, they could encounter them. To avoid such collisions during neurofeedback, the objects were rendered such that they could be passed through. Information about the state of the game on each trial, including the avatar’s location in the arena over time, was collected at 60-Hz and saved from Unity. This information was processed with custom Python scripts, downsampled to 0.5 Hz using average pooling to match the temporal resolution of the fMRI acquisitions, and shifted 2 TRs (4 s) to correct for the hemodynamic lag.Joystick task specificationsThe joystick (Current Designs Tethyx) had a two-axis (X, Y ) shaft controller with an angular range of 30◦ (±15◦) and a square range in X and Y movement planes of ±15◦. Participants could move the avatar forward at a constant speed and rotate it at any angle (−180◦ to 180◦). Upon release, an internal spring returned the shaft to vertical. The joystick was connected to the experimental computer via a fiber optic cable and interfaced directly with Unity as a standard game controller. Participants held the joystick shaft with their dominant hand and placed the base on the scanner bed or on their abdomen, securing it with their non-dominant hand. Once situated inside the scanner, participants were allowed up to two practice trials before the first scan began. During the joystick session, the time limit was 60 s in the first run and +10 s for each subsequent run as the distance between spawning and flag locations increased. Participants were not penalized if they were unable to reach the goal. The resulting number of trials per run varied (mean = 19.1 trials, range = 8–26). A trial ended upon reaching the flag, after which the game paused for 6 s (3 TRs) before the next trial.Additional MRI parametersAfter completing the functional runs in each session, we acquired two spin-echo field map volumes in opposite phase encoding directions to perform field-map correction (TR = 5 s; TE = 80 ms; 37 axial slices; voxel size = 3 mm isotropic; flip angle = 90◦; IPAT GRAPPA acceleration factor = 2). During the joystick session, we additionally collected two high-resolution anatomical images for spatial registration: a 3-D T1-weighted MPRAGE scan (TR = 2.5 s; TE = 2.9 ms; voxel size = 1 mm isotropic; flip angle = 8◦; IPAT GRAPPA acceleration factor = 2) and a 3-D T2-weighted fast spin echo scan with variable flip angle (TR = 3.2 s; TE = 565 ms; voxel size = 1 mm isotropic; IPAT GRAPPA acceleration factor = 2).Validating neural realignmentAfter observing an increase in ∆P EV along the IM and WMP components during their respective training sessions, we designed three additional analyses to further establish neural realignment as a candidate neural strategy underlying BCI learning in our task. 2 ### Page 4 First, we asked whether there was a difference in the overall variance in brain activity data before versus after learning. We calculated the total variance across all samples, in both voxel space and T- PHATE space, in the first and last runs for each session type. We then subtracted the first from last run total variance and assessed the statistical significance of the difference using bootstrap resampling (10,000 iterations) against the null hypothesis of zero change (Figure S8A). Second, we investigated the specificity of ∆P EV (Figure 4) to the components being trained, rel- ative to arbitrary other dimensions of the neural subspace. For each session type and participant, we calculated ∆P EV along each of the 20 T-PHATE manifold components (not just the one trained with neurofeedback in that session). We created a null distribution for each participant and session by randomly sampling from these 20 ∆P EV scores 10,000 times. The means of the resulting null distribu- tions did not differ from zero across participants for the IM (P = 0.619), WMP (P = 0.856), or OMP (P = 0.431) sessions, suggesting that positive ∆P EV was not observed generically across manifold com- ponents. We then computed the z-score of each participant’s true ∆P EV along the trained component in a given session (e.g., CIM during the IM session) relative to the corresponding null distribution and compared these z-scores across participants against zero to evaluate whether ∆P EV was significantly stronger for the trained component than for arbitrary components (Figure S8B). Third, we compared neural realignment to an alternative explanation for ∆P EV . The trained com- ponent was mapped continuously to the full range of turning directions in the game, so participants needed to generate a diverse range of projections to achieve flexible control of the avatar’s movement. Thus, we hypothesized that on-manifold learning would increase the variance of projections along the trained component and that this increased variance would reflect a scaling of the distribution of projec- tions (realignment). However, an increase in variance could also occur if participants learned to generate more of the larger rather than smaller projections within the same distribution of projections. Although such a “subselection” strategy would come at the expense of being able to use smaller projections for certain turning directions, we sought to evaluate these explanations quantitatively. We first ascertained that during the manifold calibration procedure, participants did create brain activity patterns which, when projected onto CIM , CW M P , and COM P , generated a range of different movement directions (Figure S13). Then, we evaluated the variance in these projections onto neuro- feedback components before and after a given component was trained (e.g., the change in variance of projections onto CIM during the IM session) (Figure S8C) and modeled whether the observed change in variance was better fit by models of realignment or subselection. Namely, based on a participant’s pro- jections onto the neurofeedback component before training (separately for each session), we simulated how the distribution would change after training as predicted by a realignment model or a subselection model, with the proportion of change controlled by a parameter P ∈ {0.2, 0.3, 0.4, 0.5}. To model subselection, we pseudo-randomly dropped a subset of points from the initial pre-training distribution (e.g., if P = 0.2 and the distribution has 100 points, 20 points would be dropped). The probability of a point being dropped was inversely proportional to its distance from 0 (i.e., lower magni- tude points were more likely to be dropped). The dropped points were then replaced by sampling from the remaining distribution. To model realignment, we first determined the best-fit probability distribution function of the original data (selected from normal, lognorm, beta, gamma, and triangular distributions) and estimated its shape parameters. Then, we generated a new distribution sampled from the best-fit function with identical shape parameters but a wider distribution scale (by increasing the scale by a factor of P :newScale = (1 + P ) ∗ originalScale). This procedure resulted in a new set of values sampled from a similar but wider distribution. We hypothesized that these two neural strategies would make different predictions about the change in projection variance (Figure S9A). We computed the change in projection variance between the simu- lated post-training distributions and the initial pre-training distribution, repeating this procedure 10,000 times for each model and proportion change level (P ). To evaluate the models, we calculated how typical the actual variance change score was under the simulated change distributions. We used Gaussian kernel density estimation to create a smooth probability density function from the simulated scores, extracted the density value of the actual score from the fitted curve, and normalized it to the maximum density of the curve. This resulted in a score between 0 and 1, where scores near 0 indicate that the actual score was highly atypical under the strategy model and scores near 1 indicate it was highly typical. We compared normalized density scores across models using a randomization test (Figure S9B). 3 ### Page 5 Intrinsic manifold stabilityIn our study, we learned the T-PHATE manifold on data from the first session (during the joystick task) and then used this manifold as the basis for the subsequent neurofeedback sessions. This approach as- sumed that the intrinsic manifold was stable across sessions. To test this assumption, we re-estimated the participant’s intrinsic manifold in each of the three neurofeedback sessions and computed the Gromov- Wasserstein distance of these manifolds from their intrinsic manifold in the first session. This distance metric is an extension of the Wasserstein distance/Earth Mover’s Distance that does not assume that the underlying spaces of two data distributions are the same. To establish stability, we assessed whether a participant’s manifold in one session was more similar to (or less distant from) their own manifold in another session than to other participants’ manifolds during the joystick or neurofeedback sessions. Let D(X1,X2) refer to the Gromov-Wasserstein distance (dgw) between two matrices, M n1 refer to participant n’s intrinsic manifold calculated from the first session, and M n i refer to the participant’s intrinsic manifold calculated from session i, which is one of their neurofeedback sessions (i.e., second, third, or fourth sessions). Thus, the within-participant stability of intrinsic manifolds is calculated as:D(n,n) = 1 34Xi=2dgw (M n1 , M n i ) To interpret this distance, we established two null distributions against which to compare D(n,n). First, we tested whether participants were significantly less distant from themselves across sessions than from others across sessions. To do this, we considered D(n,n) relative to the set of D(n,j) scores for all pairs of participants (n, j), where D(n,j) = dgwM n1 , M j i. Within each participant n, we sampled their scores D(n,j) 10,000 times to establish a null distribution and computed the z-score of D(n,n) relative to this null distribution to test whether participants were significantly less distant from themselves across sessions than from others across sessions. The z-scores were then compared across participants against zero using a nonparametric randomization test (two-tailed). In a more stringent test, we considered whether the manifolds of participants performing the joystick task were less distant from themselves in other (neurofeedback) sessions than from other participants performing the joystick task. This addresses whether stability is greater within-participant across ses- sions than within-session across participants. We defined this null distribution as the set of D(n,j) scores for all pairs of participants (n, j), where D(n,j) = dgwM n1 , M j1, sampled 10,000 times. As in the prior analysis, we computed the z-score of D(n,n) relative to this null distribution and statistically compared the resulting z-scores against zero (Figure S14). 4 ### Page 6 Supplementary ResultsParticipant Self-reported strategies5 (1) Thinking about moving straight, (2) thinking about following the yellow line, (3) following the diagonal path along the green check marks, (4) thinking about heading directly to the flagpole (most successful). 6 (1) Thinking about nice things (vacations, dates, etc), (2) doing arithmetic, (3) planning for the next months, (4) looking at the flag, (5) looking at the avatar’s feet (most successful). 7 Various strategies, none of which were obviously successful. 8 (1) Thinking about happy, sad, stressful, and pleasurable memories ((pleasur- able memories were more successful than non-pleasurable), (2) replaying movies that they had seen multiple times and could recall vivid visual and auditory features (repeating the same movie was less effective the second time), (3) thinking about the octopus on the flag (successful some of the time). 10 (1) Thinking about walking, (2) thinking about playing games, (3) thinking about playing sports, (4) imagining moving toward the flag, (4) focusing on the scenery and trees, (5) reminiscing old memories, (6) trying to embody happiness, sadness, or other emotions. 11 Playing a song in my head in Swahili. If the song was in Swahili, the avatar went straight to the flag, but if I accidentally mentally sang songs in English, the avatar spun. I also instructed the avatar at the same time, but the instructions had to be in English while the song was still in Swahili. This was very successful but if I scrambled my languages, it took a while to get working again; overall, it was easy to overthink and get confused. 13 (1) Thinking “straight” to keep the avatar on the straight path (most success- ful), (2) keeping my eyes on the goal location, (3) envisioning a sprinter running down the final stretch of a race, (4) envisioning myself walking, (5) reminding myself to not get distracted, (6) thinking “right” or “left” to get back on the path when the avatar veered off course. 14 (1) Singing to myself (not successful), (2) thinking happy thoughts (somewhat successful), (3) pretending I was holding the controller (somewhat successful), (4) pretending I was walking to a destination while imagining and creating a mental scene (most successful). 15 (1) Looking for patterns in the squares ahead, (2) imagining I was walking my dog and he was pulling me toward the flag because there was food or other dogs over there (most successful but required very vivid imagery). 16 (1) Following target with eyes and repeating “octopus flag” over and over (successful in the beginning then stopped), (2) following the target with eyes and thinking about daily life (most successful), (3) imagining moving the joy- stick with hand (sometimes successful), (4) grinding teeth when going off path (painful and not successful). 17 (1) Telling the avatar mentally to move, (2) relaxing my mind and eyes by focusing on the sound of the MRI and trying not to get stressed out by going in the wrong direction (most successful). 18 (1) Thinking “forward” over and over (somewhat successful initially at low levels), (2) thinking about my mind moving the avatar to the target, (3) staring intensely at the flag (most successful). 19 (1) Thinking “left” and “right” (most successful), (2) thinking about tapping on my leg. 5 ### Page 7 21 (1) Imagining running to break up a fight between a pair of kids before one gets bitten, (2) thinking about walking to work, (3) mentally running through my to-do list, (4) relaxing and counting the avatar’s steps, either to get to the flag or 1 to 10 on repeat (most successful). 22 Thinking about cooking different meals from start to finish. 23 (1) Imagining scenarios that involved forward, one-directional movement like dog walking, (2) imagining I was surfing or swimming (most successful). 24 (1) Pretending that someone or something I wanted to reach was at the flag, (2) imagining I was walking down the street to my house (both were successful). 25 (1) Mentally singing “Just Keep Swimming” from Finding Nemo, (2) imagining walking to class, (3) intentionally not focusing on the screen and trying not to overthink my performance (most successful). Table 1: Participant strategies. Participants (N = 18) completed a debriefing questionnaire after all sessions. They were asked: “What strategies did you try in order to control the avatar? How successful did you think these strategies were?” They typed their answers into a blank text box. Responses here are lightly edited for clarity and conciseness. 6 ### Page 8 z=-28 z=-12 z=0 z=18 z=32 z=48 z=64y=-98 y=-74 y=-60 y=-40 y=-16 y=2 y=30x=46 x=34 x=22 x=10 x=-2 x=-14 x=-26Figure S1: Navigation network targeted with neurofeedback. Sagittal, coronal, and axial slices of the MNI template brain showing voxels that were used for neurofeedback. Primary regions of in- terest include the hippocampus, medial temporal lobe cortex, precuneus, lateral occipital cortex, and parahippocampal cortex. 7 ### Page 9 PC 1PC 2T-PHATE 1 T-PHATE 2 PC 1T-PHATE 1 PC 1T-PHATE 1 PC 1T-PHATE 1PC 1T-PHATE 1B Embeddings with PCAA Embeddings with T-PHATEParticipant 5 Participant 7 Participant 13 Participant 24 Participant 25 Participant 5 Participant 7 Participant 13 Participant 24 Participant 25 LEGENDPoints: brain activity in embeddingColor: corresponding arena locationx coordinate y coordinate2 3 5 10 15 20 # embedding dimensions 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 Mean-squared error 2 3 5 10 15 20 # embedding dimensions -5 0 5 10 15 20 25 z-scoreEmbedding type T-PHATE PCAC Decoding location in arena D Arena-brain representational similarityFigure S2: Comparing T-PHATE manifolds with PCA embeddings. (A) Two-dimensional T- PHATE embeddings of the same participants in Figure 1B, showing timeseries during four consecutive trials of the joystick with timepoints colored by the coordinate location of the avatar in the game arena (color map shown in Legend). Visually, T-PHATE embeddings show clear trajectories through the latent space and clustering of nearby locations within the game arena. (B) Two-dimensional PCA projections of the exact same data. Data in PCA space are more scattered, not capturing trajectories through the latent space or much clustering of nearby locations within the game arena. These visualizations support that T-PHATE yields more meaningful embedding dimensions than PCA. (C) Mean-squared error of linear models trained to predict the location of the avatar in the game arena at each time point from data embedded in [2, 3, 5, 10, 15, 20] dimensions with either PCA or T-PHATE. T-PHATE-based decoding error is stable across multiple embedding dimensionalities and lower than PCA-based decoding error. (D) Representational similarity analysis between distances in the game arena and distances of neural embeddings with T-PHATE or PCA. T-PHATE embeddings show higher representational similarity, suggesting they better represent task dynamics than PCA embeddings, and in fewer dimensions. In (C) and (D), dots reflect the mean across participants and error bars reflect the standard error of the mean. 8 ### Page 10 B Average manifold regularization penalty0 10000 0 2x10 5 1x10 5 Participants Simulated 2x10 -21x10 -1 0 0 10000 Mean-squared error epochs epochsA Average MRAE training loss0 5x10 6 1x10 7 0 100002.5x10 -1 5x10 -1 7.5x10 -1 1x10 0 0 10000 Participants Simulated Mean-squared error epochsepochsFigure S3: Average training losses. (A) The training loss for MRAE combines a reconstruction penalty plus a manifold regularization penalty, shown for real participants (left) and simulated fMRI data as a null comparison (right). (B) Plotting the manifold regularization penalty component of the training loss separately. Loss functions are visualized as the mean (line) and 95% confidence interval (band) of the mean across participants. Individual participant curves are visualized in Figures S4 and S5. 9 ### Page 11 -0.5 0.0 log training loss avatarRT_sub_05 avatarRT_sub_06 avatarRT_sub_07 avatarRT_sub_08 avatarRT_sub_09 -0.5 0.0 log training loss avatarRT_sub_10 avatarRT_sub_11 avatarRT_sub_13 avatarRT_sub_14 avatarRT_sub_15 -0.5 0.0 log training loss avatarRT_sub_16 avatarRT_sub_17 avatarRT_sub_18 avatarRT_sub_19 avatarRT_sub_20 0 10000 -0.5 0.0 log training loss avatarRT_sub_21 0 10000 avatarRT_sub_22 0 10000 avatarRT_sub_23 0 10000 avatarRT_sub_24 0 10000 avatarRT_sub_25 2 4 6 8 log training loss sim_0 sim_1 sim_2 sim_3 sim_4 2 4 6 8 log training loss sim_5 sim_6 sim_7 sim_8 sim_9 2 4 6 8 log training loss sim_10 sim_11 sim_12 sim_13 sim_14 0 10000 2 4 6 8 log training loss sim_15 0 10000 sim_16 0 10000 sim_17 0 10000 sim_18 0 10000 sim_19 Neurofeedback participants Simulated participantsOverall MRAE training lossepoch epoch epoch epoch epoch epoch epoch epoch epoch epochFigure S4: Participant-level MRAE training loss. MRAE training loss curves (L) of 10,000 training epochs for real participants (top, black) and simulated null participants (bottom, gray). See Figure S3A for averages. Overall loss is the weighted combination of the reconstruction error and the manifold regularization penalty. MRAE training loss decreased rapidly and plateaued for all real participants by ∼ 5, 000 epochs, though final training loss plateaued at different points across participants. As expected, simulated null participants showed no decrease or convergence (note the substantially higher y-axis range). 10 ### Page 12 2 4 6 log embedding loss sim_0 sim_1 sim_2 sim_3 sim_4 2 4 6 log embedding loss sim_5 sim_6 sim_7 sim_8 sim_9 2 4 6 log embedding loss sim_10 sim_11 sim_12 sim_13 sim_14 0 10000 2 4 6 log embedding loss sim_15 0 10000 sim_16 0 10000 sim_17 0 10000 sim_18 0 10000 sim_19 Neurofeedback participants Simulated participants -3.0 -2.0 log embedding loss avatarRT_sub_05 avatarRT_sub_06 avatarRT_sub_07 avatarRT_sub_08 avatarRT_sub_09 -3.0 -2.0 log embedding loss avatarRT_sub_10 avatarRT_sub_11 avatarRT_sub_13 avatarRT_sub_14 avatarRT_sub_15 -3.0 -2.0 log embedding loss avatarRT_sub_16 avatarRT_sub_17 avatarRT_sub_18 avatarRT_sub_19 avatarRT_sub_20 0 10000 -3.0 -2.0 log embedding loss avatarRT_sub_21 0 10000 avatarRT_sub_22 0 10000 avatarRT_sub_23 0 10000 avatarRT_sub_24 0 10000 avatarRT_sub_25MRAE Manifold regularization penaltyepoch epoch epoch epoch epoch epoch epoch epoch epoch epochFigure S5: Participant-level manifold regularization penalty. Manifold regularization penalty curves (Lgeometric) of 10,000 training epochs for real participants (top, black) and simulated null partic- ipants (bottom, gray). See Figure S3B for averages. The decreasing values for real participants indicate that the encoder f learned weights that successfully embedded new samples into their respective posi- tions along the T-PHATE manifold, meaning that there was task-related and predictable structure in the input data that matched the T-PHATE representation. This was not possible for the simulated null participants, who showed no decrease or convergence (note the substantially higher y-axis). 11 ### Page 13 WMP OMP Session type -80 -40 0 40 80Δ Brain Control (%) n.s. n. s. n. s.WMP firstOMP first ** Session orderEffect of session order on BCI learning**Figure S6: Effect of perturbation order on BCI learning. Participants received the WMP and OMP in two separate sessions, with the ordering of these sessions counterbalanced across par- ticipants. Dots correspond to ∆BrainControl for individual participants by session type and sorted based on whether they completed WMP (pink) or OMP (red) first. Bars correspond to the average ∆BrainControl across participants per session type and order; error bars correspond to 95% confidence interval of the mean. We analyzed performance differences using an ordinary least squares (OLS) regres- sion with session type (within participant) and session order (between group) as factors and their inter- action term. We found a significant main effect of session type (F (1, 32) = 9.18, P = 0.005, η2p = 0.22), no main effect of session order (F (1, 32) = 1.03, P = 0.32, η2p = 0.03), and no interaction between session type and session order (F (1, 32) = 1.14, P = 0.29, η2p = 0.03). Thus, participants did not reliably differ in their ability to learn depending upon which perturbation they received first. Nevertheless, WMP learning was numerically clearer when that session came first after IM compared to when it followed OMP, perhaps consistent with a form of “learned helplessness” or a continuation of ineffective strate- gies. That said, these analyses should be interpreted with caution because testing counterbalancing order required dividing the sample size in half, reducing statistical power. ** P < 0.01, n.s. P ≥ 0.1 12 ### Page 14 IM WMP OMP IM WMP OMPIM WMP Component measured Component measured Component measuredIM Session OMP SessionWMP SessionΔ Explained variance (%)Change along all components during a specific component’s training session-10 -5 0 5 10 ** n.s. n.s.~n.s.n.s.OMn.s. * n.s.*n.s. **n.s. n.s. n.s. n.s. n.s.n.s.Figure S7: Neural realignment with all feedback components. We measured the change in the percentage of explained variance (∆P EV ) along the IM, WMP, and OMP components over the course of training for each component. That is, for the IM training session, we measured the ∆P EV for the IM, WMP, and OMP components, though the WMP and OMP components were not being trained during that session. Dots correspond to the ∆P EV for individual participants per component and session type; lines connect dots from the same participant. Bars represent the mean across participants per component and session type. Neural realignment was specific to the component being trained, as IM and WMP showed positive ∆P EV only when they were being trained in the IM and WMP sessions, respectively; as expected, none of the components showed reliable ∆P EV in the OMP session in which learning did not occur. For the IM session, mean ∆P EV along: CIM = 1.06, 95% CI = [0.35, 1.72], P = 0.005, d = 0.70;CW M P = 0.26, 95% CI = [−0.48, 0.99], P = 0.53, d = 0.15; COM P = 0.03, 95% CI = [−1.20, 1.49], P = 0.96, d = 0.01. For the WMP session, mean ∆P EV along: CIM = −0.40, 95% CI = [−1.30, 0.42], P = 0.39, d = −0.21; CW M P = 0.71, 95% CI = [0.01, 1.49], P = 0.04, d = 0.42; COM P = −1.18, 95% CI = [−2.66, 0.06], P = 0.13, d = −0.38. For the OMP session, mean ∆P EV along: CIM = −0.38, 95% CI = [−1.51, 0.92], P = 0.58, d = −0.14; CW M P = 0.36, 95% CI = [−0.77, 1.43], P = 0.55, d = 0.15; COM P =−0.45, 95% CI = [−1.36, 0.39], P = 0.35, d = −0.23. ** P < 0.01, * P < 0.05, ∼ P < 0.1, n.s. P ≥ 0.1 . 13 ### Page 15 Change in overall navigation network activity variance pre/post learningIM WMP OMP Session type -2 -1 0 1 2 3 4 Z-score** * n.s.AIM WMP OMP -0.8 -0.4 0.0 0.4 0.8Voxel resolution IM WMP OMPT-PHATE embeddingΔ total variance ~ ~ n.s.Change in component projectionsCΔ PEV vs. random componentsBn.s. n.s. n.s. n.s. n.s. n.s. Session type Session type IM WMP OMP Session type-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3Δ projection varianceFigure S8: Evaluating alternative explanations of neural variance changes. (A) For each session type and participant, we computed the total variance across all samples and voxels (left) or T-PHATE components (right) and subtracted the first run’s total variance (pre) from the last run’s total variance (post). There was no evidence of change in total variance (P s ≥ 0.17). (B) We calculated ∆P EV along each of the 20 manifold components and created a null distribution for each participant and session type by randomly sampling from these values 10,000 times. There was no evidence of positive ∆P EV for arbitrary dimensions in any session type based on the mean of these distributions evaluated across participants (P s ≥ 0.431). Next, we calculated the z-score of the participant’s actual ∆P EV for the trained component in a given session relative to the corresponding null distribution. The z-scored ∆P EV values were significantly different from zero for CIM (mean z = 0.57, 95% CI = [0.25, 0.88], P = 0.002, d = 0.81) and CW M P (mean z = 0.41, 95% CI = [0.03, 0.82], P = 0.031, d = 0.47), but not COM P (mean z = −0.17, 95% CI = [−0.59, 0.24], P = 0.77, d = −0.18). (C) We quantified how well participants could learn to generate avatar turning directions by calculating the variance of the magnitude of projections onto the trained component (e.g., along CIM during the IM session). There was a marginal increase in projection variance from the first to last run for the IM (mean difference = 0.028, 95% CI = [0.00, 0.06], P = 0.06, d = 0.38) and WMP (mean difference = 0.013, 95% CI = [0.00, 0.03], P = 0.079, d = 0.35) sessions, but not in the OMP session (mean difference = −0.027, 95% CI = [−0.06, 0.00], P = 0.95, d = −0.04). Dots correspond to individual participants; lines connect dots from the same participant. Bars represent the mean across participants. ** P < 0.01, * P < 0.05, ∼ P < 0.1, n.s. P ≥ 0.1. 14 ### Page 16 A Schematic: How strategy model and proportion of points reassigned impact projection variance B Typicality of observed Δ projection variance under each strategy modelProportion reassigned Density (normalized)IM SessionProportion reassignedWMP Session-1.0 -0.5 0.0 0.5 1.0 Projection value 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 Density Proportion reassigned = 0.2 original var=0.091 subselection var=0.130 realignment var=0.193originalsubselectionrealignment-1.0 -0.5 0.0 0.5 1.0 Projection value Proportion reassigned = 0.5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 Density original var=0.091 subselection var=0.097 realignment var=0.109 original subselection realignment subselection realignmentStrategy model Strategy modelStrategy model0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1.0 0.2 0.3 0.4 0.5Figure S9: Modeling possible neural strategies. (A) We tested whether the overall increase in variance of projections reflected either subselection (orange), where smaller projections become less frequent, or realignment (purple), where the range of projections becomes wider. We modeled these changes by taking the original distribution from the first run of a session (green) and simulating changed distributions representing different degrees of each strategy (controlled by PR). (B) (Left) During the IM session, the realignment model provided a significantly better fit to the observed scores than the subselection model at the lowest value of PR; at higher values, the difference was similar numerically but less reliable (PR = 0.2: difference = 0.21, 95% CI = [0.06, 0.37], P = 0.017, d = 0.62; PR = 0.3: difference = 0.20, 95% CI = [−0.04, 0.42], P = 0.116, d = 0.39; PR = 0.4: difference = 0.18, 95% CI = [−0.07, 0.42], P = 0.187, d = 0.33; PR = 0.5: difference = 0.19, 95% CI = [−0.06, 0.43], P = 0.165,d = 0.35). These weaker effects could relate to the finding that a participant’s ∆P EV did not predict their ∆BrainControl in the IM session: given that the IM session was always first and CIM was most geometrically consistent with the intrinsic manifold, learning to control CIM could be supported by a mixture of these strategies rather than the optimal solution of realignment. (Right) During the WMP session, the realignment model predicted the observed changes significantly better than the subselection model at all levels of PR, with the difference growing as PR increased (PR = 0.2: difference = 0.26, 95% CI = [0.03, 0.48], P = 0.038, d = 0.52; PR = 0.3: difference = 0.33, 95% CI = [0.11, 0.55],P = 0.0098, d = 0.68; PR = 0.4: difference = 0.40, 95% CI = [0.21, 0.60], P = 0.001, d = 0.94;PR = 0.5: difference = 0.42, 95% CI = [0.14, 0.62], P = 0.0006, d = 1.00). Thus, these model simulations bolster the conclusion that neural realignment enabled successful learning during the WMP session, and to a lesser extent during the IM session. Dots correspond to the normalized density for individual participants; lines connect dots from the same participant. Bars represent the mean across participants. *** P < 0.001, ** P < 0.01, * P < 0.05, ∼ P < 0.1. 15 ### Page 17 Repeat for first and final runsVoxels ...β0 β1 ... β6 β0 β1 ... β6XTimepoints Training set K-1 trialsTrial K Testing set apply weightsβ0 β1 ... β6 β0 β1 ... β6X = = x coordinate y coordinate Trials 1:K-1 coordinates Trial K coordinates x coordinate y coordinateprediction error = distance(predicted, true) predicted true Repeat for K folds Center searchlight on each voxel in the brainSearchlight-based location decodingFigure S10: Whole-brain searchlight analysis. An exploratory searchlight analysis was used to test for improved task representation in regions across the whole brain, by decoding the avatar’s coordinates in the game arena from each timepoint using ridge regression with nested cross-validation. Regression models were trained and tested with leave-one-trial-out cross-validation in the outer loop, and hyper- parameter (α) optimization was performed in the inner loop (α ∈ {10k : k ∈ {−5, −4, −3, . . . , 5}}). Regression models were scored as the mean squared error (MSE) between model-predicted arena loca- tions and true participant locations. We measured decoding improvement from the voxelwise map of MSE across the brain in the first minus last runs of a session; lower MSE means improved prediction and thus a clearer representation of the task. 16 ### Page 18 JI HGFPC 1 PC 2 Diffusion component 1 Diffusion component 2 PC 1 PC 2 Diffusion component 1 Diffusion component 2 T-PHATE 1 T-PHATE 2 PCA embedding Diffusion map embedding T-PHATE embedding T-PHATE 2 T-PHATE 1 T-PHATE embeddingDiffusion map embeddingPCA embeddingA ELKGround truth labelsGround truth labels3D Swiss roll, no noise 3D Swiss roll, high noisedim 1 dim 2Ground truth labels dim 2 dim 1 dim 3 2D Swiss roll, high noise dim 2 dim 1 dim 3 dim 1 dim 2PC 1 dim 1 dim 2Diffusion component 1 dim 1 dim 2T-PHATE Eig 1 Labeled by PC 1 Labeled by diffusion component 1 Labeled by T-PHATE Eigenvector 1B C DFigure S11: Nonlinear methods improve detection of on-manifold activity in nonlinear data. (A) Simulated, noise-free Swiss roll dataset of 1,000 points displayed in ambient 3-D space (standard deviation of Gaussian noise, SD = 0.0). Data colors represent position along the main dimension of the ground-truth manifold. (B–D) 2-D embeddings of (A) with (B) PCA, (C) diffusion maps, and (D) T-PHATE. Nonlinear dimensions (i.e., diffusion maps and T-PHATE) reveal an “unrolled” representation of the Swiss roll data, representative of its structure. (E) Swiss roll dataset with high noise (SD = 1.0). (F–H) Embeddings of (E) as in (B–D). With higher noise added to the 3-D data, T-PHATE remains able to unroll the true data manifold, whereas the diffusion map remains partially coiled. (I) Swiss roll with high noise in a 2-D ambient space, colored by its ground-truth manifold positions. (J–L) Data as in (I), now colored by each point’s loading onto the first component of the embeddings from (J) PCA, (K) diffusion maps, and (L) T-PHATE. As in the ground truth labeling (I), the loadings shown in (L) track the coil of the Swiss roll, whereas (J, K) do not.17 ### Page 19 0 10 20 30 40 50 % Variance Explained avatarRT_sub_05 avatarRT_sub_06 avatarRT_sub_07 avatarRT_sub_08 avatarRT_sub_10 avatarRT_sub_11 0 10 20 30 40 50 % Variance Explained avatarRT_sub_13 avatarRT_sub_14 avatarRT_sub_15 avatarRT_sub_16 avatarRT_sub_17 avatarRT_sub_18 5 10 15 20 Component 0 10 20 30 40 50 % Variance Explained avatarRT_sub_19 5 10 15 20 Component avatarRT_sub_21 5 10 15 20 Component avatarRT_sub_22 5 10 15 20 Component avatarRT_sub_23 5 10 15 20 Component avatarRT_sub_24 5 10 15 20 Component avatarRT_sub_25 Eigenspectra of participants’ intrinsic manifolds X = IM, WMP, or OMP = untrainedFigure S12: Evaluating eigenspectra of participants’ intrinsic manifolds. For each participant, we calculated the percentage of variance explained by each of the 20 eigenvectors of their intrinsic manifold. Components that were trained (1st, 2nd, and 20th components for each participant; CIM ,CW M P , and COM P , respectively) are shown as red X markers; untrained components are visualized as blue points. The difference in percentage of variance explained by CIM and CW M P varied across individual participants, but we did not find any clear relationships between a participant’s difference in percentage of explained variance and neural (∆P EV ) or behavioral (∆BrainControl) learning. 18 ### Page 20 0.0 0.5 1.0 1.5 2.0 Density avatarRT_sub_05 avatarRT_sub_06 avatarRT_sub_07 avatarRT_sub_08 avatarRT_sub_090.0 0.5 1.0 1.5 2.0 Density avatarRT_sub_10 avatarRT_sub_11 avatarRT_sub_13 avatarRT_sub_14 avatarRT_sub_150.0 0.5 1.0 1.5 2.0 Density avatarRT_sub_16 avatarRT_sub_17 avatarRT_sub_18 avatarRT_sub_19 avatarRT_sub_20-90 0 90 0.0 0.5 1.0 1.5 2.0 Density avatarRT_sub_21-90 0 90 avatarRT_sub_22-90 0 90 avatarRT_sub_23-90 0 90 avatarRT_sub_24-90 0 90 avatarRT_sub_25 Decoded angle Decoded angle Decoded angle Decoded angle Decoded angle IM WMP OMPDistributions of joystick session testing data projected onto feedback componentsFigure S13: Projections onto neurofeedback components. With the 20% of testing data from the joystick task used to validate our decoding pipeline, we computed the loadings of the brain activity data onto the MRAE subspace, projected them onto the three feedback components, and then mapped those values to an avatar turning direction between −90◦ and 90◦. Distributions of these directions (x-axis) are shown separately for each participant, with the same test data from their first session projected onto the IM (navy), WMP (teal) and OMP (green) components in their T-PHATE manifold. 19 ### Page 21 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 z-score distance *** -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 -2.5 -2.5 ***Within-participant intrinsic manifold stability across sessionsSelf vs. others, across sessions Self vs. others, within session 1Figure S14: Evaluating stability of intrinsic manifolds across sessions. The stability of the T- PHATE manifold across sessions was assessed using Gromov-Wasserstein distance. The distance between a participant’s manifold learned from the joystick session and in the neurofeedback sessions (D(N,N ) was compared against: (left) the distances between the participant’s joystick manifold and the neurofeedback manifolds of other participants (mean z = 1.051, P = 0.0002, 95% CI = [−1.3422, −0.7168], d = −1.49), and (right) the distances between the participant’s joystick manifold and the joystick manifold of other participants (mean z = 0.416, P = 0.0006, 95% CI = [−0.55, −0.27], d = −1.34). Thus, manifolds were significantly more similar within vs. between participants, confirming their stability over sessions. Dots correspond to the z-score of each individual’s within-participant distance relative to the distribution of between-participant distances. Bars represent the mean across participants; error bars reflect the 95% CI of the mean. *** P < 0.001 20 ### Page 22 WMP OMP Session type -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 n.s.~IM WMP OMP Session type 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 MSE n.s. n.s. n.s. Decoders trained/tested within session Δ error of IM-trained decodersA BΔ MSEFigure S15: Evaluating decoders across sessions. Decoding of game arena location was performed as in Figure 1C. (A) Decoders trained and tested within session with leave-one-run-out cross-validation performed equally well in terms of mean squared error (MSE) (P s ≥ 0.2), suggesting that brain activity from all three sessions equally represented the arena (see Figure S2C for comparison with voxel-space and PCA decoding). (B) Decoders were trained on all runs from the IM session and then tested on the first and last runs of the WMP and OMP sessions, respectively. Change in decoding accuracy was calculated by subtracting the MSE for the first run from the MSE of the last run (higher reflects more error after training). IM-trained decoders performed marginally worse after WMP training (mean ∆M SE = 0.049, P = 0.0501, 95% CI = [−0.001, 0.106], d = 0.41) but not after OMP training (mean ∆M SE = 0.007, P = 0.369, 95% CI = [−0.0370, 0.048], d = 0.076). Combined, these results suggest that the arena was represented similarly in brain activity during the IM and WMP sessions, but the subspace defined by WMP diverged from IM after successful learning; this divergence was not observed in the OMP condition, perhaps because learning was unsuccessful. Dots correspond to individual participants. Bars represent the mean across participants; error bars reflect the 95% CI of the mean. ∼ P < 0.1, n.s.P ≥ 0.1. 21 ### Page 23 B Correlations between mask size and learning metrics 05 06 07 08 10 11 13 14 15 16 17 18 19 21 22 23 24 25 Participant ID 0 Count 400 800 1200 1600AIM WMP OMP-1.00 -0.50 0.00 0.50 1.00 BCI learning (ΔBrainControl) IM WMP OMP Neural learning (ΔPEV)Spearman’s ρn.s. n.s.n.s. n.s. n.s. n.s. Number of voxels included in navigation maskFigure S16: No relationship between size of neurofeedback mask and successful learning.(A) After aligning the navigation mask from Neurosynth.org into native space, participants varied in the number of voxels retained for neurofeedback (mean = 1, 354 voxels; range = 1, 058 − 1, 599). (B) (Left) No significant relationship was found between a participant’s mask size and their behav- ioral learning (∆BrainControl) for any session type (Spearman’s ρ, 95% CI for each session: IM = −0.29, [−0.71, 0.22]; WMP = −0.27, [−0.72, 0.25]; OMP = 0.035, [−0.46, 0.48]; P s ≥ 0.32). (Right) No significant relationship was found for neural learning (∆P EV ) either (IM = 0.06, [−0.42, 0.53]; WMP = 0.15, [−0.33, 0.61]; OMP = 0.12, [−0.45, 0.63]; P s ≥ 0.21). Bars reflect ρ for each session; error bars represent the 95% CI of ρ from bootstrap resampling (10,000 iterations). 22