Brain Harmony: A Multimodal Foundation Model Unifying Morphology and Function into 1D Tokens

Distinction from prior foundation models

We appreciate the reviewer’s feedback and would like to clarify this important point. To the best of our knowledge, we are indeed the first to deeply compress both brain morphology and dynamics into unified representations via 1D tokens as a multimodal brain foundation model. Regarding “prior foundation models”, we would greatly appreciate it if the reviewer could point out specific works we might have overlooked.

As acknowledged by the other three reviewers, our proposed BrainHarmonix is a significant conceptual step forward, addressing a clear need in neuroimaging by effectively handling heterogeneous TRs and multimodal data fusion, thus providing novel directions for future research in brain foundational models.

Moreover, Table 1 and 2 in our paper provide how the proposed model differs from previous brain foundation models. For further clarity, we summarize the key distinctions below:

We hope this comparison clarifies the unique contributions of our work.

Data diversity

We thank the reviewer’s insightful comment and would like to emphasize that we have explicitly acknowledged the limitations regarding age coverage in our pretraining datasets within the "limitations and future work" section.

Importantly, our current pretraining data already extends the lifespan coverage compared to prior state-of-the-art brain dynamics models such as BrainLM and Brain-JEPA, expanding from middle-aged and older adults to children as well, thereby enabling downstream analyses of neurodevelopmental disorders that previous models did not address. Moreover, our pretraining datasets inherently include scanner diversity, for example, ABCD data was acquired from multiple vendors (Siemens, GE, Philips). To further demonstrate generalizability across different geographical and demographic contexts, we have extended our downstream evaluations to include an Asian clinical cohort collected from a university hospital memory clinic, effectively assessing our model in real-world clinical scenarios. Please refer to our response to Reviewer Eo3b, "Limited generalization exploration beyond public Western datasets" for the detailed results.

Evaluation on additional datasets and tasks

We thank the reviewer for this valuable suggestion. Due to the limited rebuttal time and resources, we apologize that we are unable to complete the downloading, preprocessing, and evaluation of our model on the suggested datasets. Nevertheless, we extended our evaluation to an Asian clinical cohort collected from a university hospital memory clinic, thereby assessing generalizability to non-Western populations and in real-world clinical scenarios, as correctly pointed out by the reviewer.

Specifically, we performed an additional task not included in the current version - classification of amyloid-positive/negative cognitively normal participants, which holds significant clinical value for AD prognosis and early intervention. Please refer to our response to Reviewer Eo3b, "Limited generalization exploration beyond public Western datasets" for detailed results. We will incorporate the suggested HCP-EP and TCP for downstream evaluation in our future work.

Fine-tuning evaluation across data portions

We appreciate the reviewer’s suggestion regarding few-shot evaluation. To address this, we conducted additional analyses by scaling the fine-tuning dataset using increasing proportions (20%, 40%, 60%, 80%, and 100%). The results regarding accuracy (%) are shown in the table below. Our results demonstrate a clear and consistent scaling of performance with increasing data portions. Notably, compared with prior leading baseline BrainMass 59.35% on AbideII and 65.99% on ADHD, BrainHarmonix achieves state-of-the-art performance even when fine-tuned on only 80% of the dataset, highlighting the efficiency and effectiveness of our pretrained representations. We will include this analysis in the revised version.

Scaling behavior with increasing token numbers

We appreciate your valuable comments. In our preliminary experiments, we observed that the performance tends to stabilize when the number of tokens reaches around 256. For completeness, we have included results with 512 and 1024 tokens as references. As shown in the table below, the accuracy (%) remains relatively stable beyond 256 tokens, confirming our initial observation.

Scalability and efficiency evaluation

We thank the reviewer for the valuable suggestions. We have added scalability experiments by varying the size of the harmonizer model from 22M parameters to 307M parameters. The results regarding accuracy (%) in the table below show performance improvement from 22M to 86M, while the gain from 86M to 307M is only marginal.

On the other hand, we investigated the effect of using different portions of the pretraining dataset. Specifically, we applied identical sampling proportions to both the UKB and ABCD datasets for pretraining. The corresponding results regarding accuracy (%) are reported in the table below. We observe that the model’s performance improves as the portion of the pretraining dataset increases.

Following the reviewer’s suggestion, we have also included the pretraining time (on 8 NVIDIA H100 GPUs (80GB)), as well as finetuning time (on 1 H100 GPU) and inference time (on 1 H100 GPU) on AbideII, corresponding to different model sizes (token counts) to provide a more comprehensive view of the computational cost in the table below. Larger model or more token counts lead to longer computing time.

Overall, our further analysis highlights a trade-off: larger models or increased token counts generally enhance performance but result in longer computation times.

Failure mode and limitation analysis

We appreciate the reviewer’s suggestion to analyze potential failure modes and areas of underperformance. In our current evaluations, one notable case where BrainHarmonix underperforms is on the ADHD-200. Its F1 is slightly lower than BrainHarmonix-F. This is likely due to motion artifacts in T1, as ADHD patients exhibit increased head motion during MRI acquisition. Such motion introduces noise and negatively affects structural data quality, potentially reducing multimodal fusion performance. Future work will explore methods to improve robustness against data-quality issues. On the other hand, although we have demonstrated data scaling effects in the above discussion, model performance under low-sample and few-shot learning scenarios remains an area for improvement. Future studies may address few-shot adaptation through approaches such as parameter-efficient fine-tuning or prompt-based tuning.

Lack of Baseline Against Dynamic Time Warping

We appreciate the reviewer’s valuable suggestion. However, we believe that DTW is not a suitable baseline in this setting. DTW assumes a meaningful temporal correspondence between sequences. However, in the context of resting-state fMRI, there is no ground truth temporal alignment across individuals, as each subject’s brain dynamics evolve independently and asynchronously. Therefore, applying DTW across different scans would impose artificial temporal correspondences not supported by the data.

We thank the reviewer for bringing this method to our attention. We will add the above discussion in the related work and TAPE section.

Interpretability of 1D tokens and attention maps

We appreciate the reviewer’s insightful comments regarding the interpretability of our 1D tokens. Although NeurIPS 2025 policy prohibits us from including new figures at this rebuttal stage, we would like to briefly summarize key insights from our new analyses.

We examined the attention patterns between the 128 learned 1D tokens and the modality-specific tokens (400 fMRI ROI + 1200 T1 tokens) in ASD diagnosis using ABIDE-II data. For the 400 fMRI ROI tokens, each is obtained by averaging all tokens within the corresponding ROI. We found differentiation in modality attention among the 1D tokens: 93/128 tokens attended exclusively to fMRI, 30/128 exclusively to T1, and 5/128 tokens exhibited cross-modal attention. For cross-modal tokens, we found that they exhibited key structure-function coupling such as medial prefrontal cortex in brain morphometry and default model network in brain dynamics, which have previously been demonstrated in the literature to be associated with ASD.

For the 93 fMRI-specific tokens, further analysis revealed network-level functional differentiation relevant to ASD behavioral traits. Specifically, 60/93 were network-specific, predominantly focusing on a single brain network (with >70% salient ROIs within one network), while the remaining 33 were identified as “bridge” tokens capturing interactions across multiple networks. Among the most salient network-specific tokens, temporoparietal network (implicated in social perception and language processing deficits), somatomotor network (associated with sensorimotor integration impairments), and default mode network (linked to mentalizing deficits) emerged prominently. The identified “bridge” tokens primarily captured interactions involving default, limbic, and control networks, reflecting impaired integration across sensorimotor, socioemotional, and higher-order cognitive processes - a mechanism implicated in pathophysiology of ASD. Please refer to our response to Reviewer wvBG, “Latent space analysis reflecting structural constraints” for the interpretation of the geometry-constrained fMRI latent space.

We will incorporate corresponding visualization into our revised manuscript.

Limited generalization exploration beyond public Western datasets

We thank the reviewer for raising this important point regarding generalization. To address this, we extended our evaluation to an Asian clinical cohort collected from a university hospital memory clinic (single site), thereby assessing generalizability to non-Western populations and in real-world clinical scenarios. Specifically, we performed an additional task not included in the current version - classification of amyloid-positive/negative cognitively normal participants, which holds significant clinical value for AD prognosis and intervention. As shown in the table below, BrainHarmonix achieved state-of-the-art performance in this clinically relevant, in-house setting, underscoring its robustness and cross-population generalizability. We will incorporate these additional results in the revised version.

Population-level geometric harmonics

We thank the reviewer for highlighting this important consideration. Using geometric harmonics derived from a group-level cortical surface provides robust, stable alignment by capturing common cortical geometry across individuals, thereby facilitating generalization and comparability across diverse datasets. In contrast, subject-specific cortical surfaces can be highly variable and prone to errors. However, we agree that combining group-level and individualized geometric harmonics may further enhance sensitivity and robustness, especially for atypical populations (e.g. patients) or the developing/aging brain. We will explore individualized (or age-specific) geometric harmonics in future studies which could further enhance model performance for clinical applications.

Zero-shot cross-modal transfer with harmonizer

We thank the reviewer for this insightful suggestion. While the current harmonizer requires paired multimodal inputs, a promising future direction would be to adapt it for zero-shot transfer and cross-modal imputation. One possible solution is to train modality-specific adapters, which could learn mappings from single-modality inputs to multimodal latent representations, enabling the harmonizer to impute missing modalities or perform cross-modal translation in scenarios where only one modality is available. We will add this in future work.

Synthetic testing of TRs

We appreciate this insightful suggestion. To further evaluate TAPE’s effectiveness, we tested on both the original HCP-A test set and version with samples randomly downsampled by factors of 1and 2 (equal probability, leading to TR values 1.6, 2.4). The comparable performance across conditions demonstrates TAPE's robustness.

Latent space analysis reflecting structural constraints

We thank the reviewer for suggesting this valuable qualitative analysis. As recommended, we extracted fMRI embeddings from BrainHarmonix-F and Brain-JEPA (400 ROIs, each represented by a 768-dimensional embedding) and applied t-SNE to project these embeddings onto a 2D plane. Since NeurIPS 2025 policy prohibits us from including new figures at this rebuttal stage, we performed a quantitative analysis comparing our model to Brain-JEPA. Specifically, we correlated each dimension of the t-SNE embedding with each of the 200 geometric harmonic modes across 400 ROIs. Compared to Brain-JEPA, our geometry-constrained embeddings exhibit a greater number of significantly correlated modes (p<0.05), with higher correlation strengths and significance levels on the top 5 most significant modes (see table below). On the other hand, we applied the Fisher r-to-z transformation to all correlations from the 200 harmonics for each model and conducted a two-sample t-test. Results demonstrate that correlations from our model are significantly higher overall. The correlation strength and statistical significance from the top modes, along with the overall comparison, confirm that our model is constrained by structural information more than Brain-JEPA. We will incorporate the corresponding t-SNE visualization in the revised version.

Clarification of TAPE

We thank the reviewer for this insightful question, and we appreciate the opportunity to clarify our TAPE implementation:

The matrix $B^{k^{\ast}}_k$ is defined for each TR. Specifically, for each TR ($s$) with corresponding patch size $k=round(\tau/s)$, we construct $B^{k^{\ast}}_k$ as a deterministic bilinear interpolation matrix that resizes embedding weights $\omega^{\ast}$ from the base size $k^{\ast}$ to $k$. This matrix is not learned but rather computed analytically using standard bilinear interpolation. The learnable parameters here are $\omega^{\ast}$, which are optimized via backpropagation during pretraining. The $B^{k^{\ast}}_k$ matrix provides a deterministic transformation to adapt $\omega^*$ to different patch sizes corresponding to different TR. Implementation details can also be found in the code provided in the supplementary material.

The consistency of temporal duration is guaranteed both theoretically and empirically.

Theoretical guarantee: $\omega$ from the transformation of $\omega^*$ with the pseudo-inverse of $B^{k^{\ast}}_k$ matrix $(\omega = ({B^{k^{\ast}}_k}^T)^+\omega^{\ast})$ is the minimum-norm least-square solution that minimises $E_x[(<x, \omega^{\ast}> - <B^{k^{\ast}}_kx, \omega>)^2]$, where $<\cdot,\cdot>$ denotes inner product. A full derivation could be found in PI-resize work [1].

Empirical evidence: To further evaluate TAPE’s effectiveness, we tested on both the original HCP-A test set and version with samples randomly downsampled by factors of 1 and 2 (equal probability, leading to TR values 1.6, 2.4). The comparable performance across conditions demonstrates TAPE's robustness.

Discussion on end-to-end training

We thank the reviewer for the insightful comments regarding end-to-end training. The harmonizer is already pretrained with long input sequences (i.e., 7200 fMRI tokens + 1200 T1 tokens + the number of 1D tokens), resulting in considerable memory requirements during training. Enabling end-to-end training would substantially increase GPU memory consumption and training time, potentially taking more than two weeks per run under our current infrastructure, limiting our capacity of model development and comprehensive evaluation.

However, we agree with the reviewer that jointly optimizing the unimodal encoders and the fusion module could potentially lead to further performance gains. Exploring efficient training strategies for such long-sequence transformer models, particularly in the context of neuroimaging, is a promising direction. We will discuss this important future work in the revised version.

The slightly lower F1 of BrainHarmonix compared to BrainHarmonix-F on ADHD-200 likely stems from motion artifacts in the T1 scans, especially prevalent among ADHD patients who tend to exhibit more head motion during MRI acquisition. Even though we have applied some quality control of T1, motion artifacts is unavoidable in T1 data of ADHD kids. This motion negatively impacts structural data quality, introducing noise that affects multimodal fusion. In contrast, fMRI data has gone through motion correction, of which the impact is less. We will investigate data-quality factors more systematically in future work.

Performance improvement: parameter count vs. multimodal integration

We thank the reviewer for the insightful question. We have added additional comparisons to better illustrate the performance gains from incorporating multiple modalities (2nd & 3rd columns in the table below containing the results regarding accuracy (%)) and from introducing the harmonizer module for fusion (4th-6th columns). Specifically, we concatenated the embeddings from the frozen T1 and fMRI encoders and passed them through a trainable linear layer for the classification task. Furthermore, we conducted experiments using harmonizers of different sizes, ranging from 22M parameters to 307M parameters (4th-6th columns). The results in the table below clearly demonstrate the performance improvements achieved both by adding modalities and by scaling the harmonizer module.

References:

[1] FlexiViT: One Model for All Patch Sizes. CVPR 2023.

Clarification on TAPE

We appreciate the opportunity to clarify Section 3.1.2 (lines 155-165).

In Figure 4, we illustrate the intuition behind TAPE. Given two fMRI time series with different TR ($s$), generating patches of a fixed temporal duration ($\tau$) results in varying patch sizes $k$ (where $k = round(\tau/s)$). Thus, a fixed-size embedding kernel ($\omega$) is insufficient. Lines 155-165 describe how we transform a kernel of size $k^{\ast}$ ($\omega^{\ast}$) to a kernel of an arbitrary length $k$ ($\omega$) using the pseudo-inverse of the bilinear transformation matrix $B^{k^{\ast}}_k$, as defined in Equation (2). $B^{k^{\ast}}_k$ matrix is computed analytically using standard bilinear interpolation. The learnable parameters here are $\omega^{\ast}$, which are optimized via backpropagation during pretraining. The $B^{k^{\ast}}_k$ matrix provides a deterministic transformation to adapt $\omega^*$ to different patch sizes corresponding to different TR. Implementation details can also be found in the code provided in the supplementary material.

To further improve the clarity of lines 155-165, we will first quote Figure 4 to show the intuition, and add the definition of every symbol on first use at the beginning: “Let $\tau$ be the temporal duration each patch should cover; $s$ the TR; $k$ an arbitrary patch length and $k^{\ast}$ the base patch length”. Please refer to our response to Reviewer wvBG, "Clarification of TAPE" for the theoretical and empirical evidence supporting our TAPE method.

We appreciate the reviewer’s suggestion regarding clarity improvements. Besides Section 3.1.2, please let us know if any other section require further clarifications..

Implementation and discussion of geometric harmonics

We appreciate the reviewer’s insightful questions and pointer to the recent discussion on [1].

We appreciate the reviewer’s observation regarding our positional embedding strategy. While Brain-JEPA utilizes gradient-based positional embeddings, our approach extends this by integrating both brain geometric harmonics and gradients in a complementary manner. Specifically, as detailed in lines 224–226, we average both brain geometric harmonics and gradients after linear projection for positioning. This combination enables our model to capture both multi-scale oscillatory modes and large-scale gradient axes, promoting structure-function coherence in the latent space. It leads to better downstream performance compared to gradient only positioning in Brain-JEPA as evidenced in Figure 6. We will make these implementation differences clearer in the revision.

The critiques to [1] focus on the paper’s claim that geometric harmonics, by themselves, can serve as a “winner-take-all” solution for brain dynamics reconstruction, thereby diminishing the role of the structural connectome. However, the critiques do not affect the validity of our geometric pre-alignment. The harmonics in our work are only used to provide geometry-aware positioning, we make no claim that they can fully explain/reconstruct brain dynamics. On the other hand, the harmonics are averaged with large-scale functional gradients, so it does not have a winner-take-all basis. We will acknowledge these points in our revised version. Future work can explore how structural connectome and other biological principles can be encoded into the model and whether they can further improve brain representation learning and generalizability.

Comparison to the arxiv work

We appreciate the reviewer bringing this work [2] to our attention. However, the absence of publicly available code of this work currently prevents direct comparison. It proposes a single-modality foundation model focusing exclusively on brain dynamics, similar to Brain-JEPA and BrainLM, whereas our work unifies both structural morphology and functional dynamics into compact 1D tokens. We will acknowledge this recent work in the revised version.

Grammer revision

We thank the reviewer for pointing this out. We will revise the sentence accordingly.

Evaluation of Brain-JEPA

We appreciate the reviewer for observing this. The experiments were conducted using the official Brain-JEPA implementation for both pretraining and fine-tuning, during which the Brain-JEPA encoder was fully fine-tuned. The reported results were averaged over three random data splits. This is a more robust evaluation compared to the original results in the Brain-JEPA paper, which used one fixed data split. Therefore, the performance gap could come from the splitting difference. The three data split settings were applied consistently in all of our comparisons. We hope this clarifies the observed differences in performance and will explicitly state this in the revised version.

Ablation of ABCD dataset

We thank the reviewer for this valuable suggestion. It is important to highlight that previous methods, such as Brain-JEPA, cannot readily integrate datasets with heterogeneous TR values. In contrast, our proposed TAPE module allows datasets with heterogeneous TR, allowing us to add both UKB and ABCD into the pre-training to cover both the developing and aging brain. Following the reviewer’s recommendation, we conducted an ablation on ADNI for MCI classification by pretraining BrainHarmonix-F without ABCD data (see table below). We found it still outperformed the original Brain-JEPA based on the same UKB dataset. We would also like to emphasize that the proposed data augmentation and pre-alignment techniques significantly contribute to the observed performance improvements as shown in Figure 6. We will include the ablation results in the revised version.

Clarification on subsampling times

We appreciate the reviewer’s question. Thanks to our proposed TAPE, our model uniquely handles heterogeneous TR, enabling us to introduce a novel data augmentation strategy for fMRI time series by downsampling. This approach enriches the dataset with a broader range of TR values and enlarges the sample size, significantly boosting model performance as demonstrated in Figure 6. Specifically, we downsample UKB data by factors ranging from 1 to 3, covering a TR range from ~0.7s to ~2.9s, and ABCD data by factors of 1 to 2, achieving a TR range up to 2.4s. Since single-band fMRI typically employs TR values between 2-3s and multi-band fMRI typically below 1s, further downsampling ABCD data would result in TR values exceeding 3s, which is uncommon in practical fMRI acquisitions and could lead to lower efficiency (longer training time).

Frozen unimodal encoders

We thank the reviewer for the insightful question. We chose to freeze the unimodal encoders during the fusion stage to allow the fusion module to flexibly integrate various pretrained unimodal representations. This design choice enables greater modularity and scalability, making it easier to adapt to different unimodal encoders trained under diverse settings without requiring full end-to-end retraining. As such, end-to-end training is not the main focus of this work. We would also like to emphasize that our current fusion module, even when trained with frozen unimodal encoders, already achieves significantly better performance than each unimodal model alone, demonstrating the effectiveness of our fusion strategy.

Moreover, the harmonizer is already pretrained with long input sequences (i.e., 7200 fMRI tokens + 1200 T1 tokens + the number of 1D tokens), resulting in considerable memory requirements during training. Enabling end-to-end training would substantially increase GPU memory consumption and training time, potentially taking much longer time per run under our current infrastructure, limiting our capacity of model development and comprehensive evaluation.

However, we agree with the reviewer that jointly optimizing the unimodal encoders and the fusion module could potentially lead to further performance gains. Exploring efficient training strategies for such long-sequence transformer models, particularly in the context of neuroimaging, is a promising direction and we will discuss it in future work in the revised version.

More ablation studies

We appreciate the reviewer’s thoughtful question. We initially chose ABIDE and ADHD-200 for our ablation studies because both datasets contain heterogeneous TR values. Following the reviewer’s suggestion, we additionally performed ablation studies on ADNI (see table below with results regarding accuracy (%)), where we observed a similar trend and performance pattern, reinforcing the effectiveness of our proposed model design. We will include this additional ablation in the revised version.

Figure styles

We appreciate the reviewer’s observation regarding the style of our figures. Each figure is designed specifically to best convey the particular type of information presented. Figure 1 serves as a conceptual illustration, whereas Figure 3 visualizes the geometric harmonics generated from actual data. For Figure 5, we use circle sizes to clearly indicate the number of trainable parameters, a variable that is not applicable in Figure 6. However, we agree with the reviewer that a more consistent visual design would benefit readability, and we will refine the figures accordingly in the revised version.

References:

[1] Geometric constraints on human brain function. Nature 2023.

[2] A Foundational Brain Dynamics Model via Stochastic Optimal Control. Arxiv 2025.