Let Brain Rhythm Shape Machine Intelligence for Connecting Dots on Graphs

We sincerely thank the reviewer for recognizing the novelty and importance of our work. We also greatly appreciate the constructive feedback regarding empirical performance, neuroscientific interpretability, and theoretical grounding of the proposed model. Below, we provide detailed responses to each of the reviewer’s concerns.

1. Lacking a comprehensive comparison with the latest state-of-the-art graph learning models [1-3].

[Comparison with additional graph learning methods]
We have added two baseline models as reviewer suggested: GRAND, GTN. In addition, we also included KuramotoGNN (KGNN)[1] and GraphCON [2], which are both recent graph models for oscillatory dynamics. While we initially considered GNNCert, it focuses on adversarial robustness and its reported clean accuracy is based on existing backbones GCN, GAT, and GIN, which are already part of our benchmark. Nonetheless, we would include GNNCert to the related work to acknowledge its relevance in the broader graph learning. The updated results are presented in Table below, the results show that our models still achieve significant improvements in brain datasets and retain strong performance in general graph benchmarks.

Brain states identification:

Node classification:

Graph classification:

2. The interpretability (biological plausibility) analysis is insufficient.

[More detailed interpretability (biological plausibility) analysis] We now elaborate how the discovered patterns in Figure 3 support neurobiological plausibility:

• Alignment with functional subnetworks (Figure 3a):
  BRICK captures task-specific synchronization across functionally defined brain network in HCP-A. For instance, Visual network regions (red) cluster clearly during the VISMOTOR task in the deeper representation $X^{(L)}$, indicating strong within-network coordination. These results suggest that our model respects functional boundaries in a biologically meaningful manner.
• Intersubject consistency and task-relevant encoding (Figure 3b):
  This figure shows that individuals performing the same task (e.g., “SOCIAL” in HCP-YA) form well-separated synchronization in deeper layers $X^{(L)}$. This indicates BRICK captures latent dynamics consistent across subjects, supporting the notion that neural phase alignment underlines cognitive state representation.
• Unsupervised Functional Differentiation (Figure 3c): Here, we assess whether BRICK can produce meaningful cortical parcellations without supervision. Compared to spectral clustering, BRICK yields phase-based partitions that are more spatially localized and better aligned with known functional networks [3]. This supports the claim that BRICK captures biologically coherent dynamics.

In this way, BRICK offers a mechanistic perspective grounded in neural oscillation theory, which allows us to move beyond black-box representations toward interpretable latent dynamics.

3. The paper lacks theoretical proof for its core mechanism (e.g., how oscillatory synchronisation mitigates over-smoothing). Although E is claimed to be a "natural" Lyapunov function, no detailed mathematical derivations or proofs are given.

[Why E is a natural Lyapunov function?]

1. Physical intuition

The dynamical system (Eq.(4)) $\frac{d\hat{x}\_i}{dt} = \omega_i + \gamma \phi\_i( y_i + \sum\_{j=1}^N w\_{ij} \hat{x}\_j)$ is a network of weakly-coupled oscillators driven by an external signal $y_i$. This suggests the following energy function (Eq.(5)):

$E(\hat{x}) = -\sum\_{i,j} \hat{x}\_i^\top w\_{ij} \hat{x}\_j - \sum\_i y_i^\top \hat{x}\_i$,

where the first sum represents the coupling term and the second sum is the driving term. Because the update $\frac{d\hat{x}_i}{dt}$ explicitly attempts to align $\hat{x}_i$ with its neighbors $(w\_{ij} \hat{x}\_j)$ and with $y_i$, the quantity measuring how well aligned they are, namely $E$, is the physically most natural Lyapunov candidate.

2. Lyapunov property [4]

Let $\hat{\mathbf{x}} = [\hat{x}\_1, ..., \hat{x}\_N] \in \mathbb{R}^{N \times d}$ and recall dynamical system Eq.(4) and energy $E$ Eq.(5). Because $W = [w\_{ij}]$ is symmetric (undirected),

$\nabla_{\hat{\mathbf{x}}} E = - (W \hat{\mathbf{x}} + \mathbf{y}), \quad \mathbf{y} = [y_1 \ldots y_N]$

Define $u_i := y_i + \sum_j w\_{ij} \hat{x}\_j$. Then Eq.(4) may be rewritten as:

$\dot{\hat{x}}_i = \omega_i + \gamma \phi\_i(u_i) = \omega\_i - \gamma \phi\_i(-u\_i) = \omega\_i - \gamma \phi\_i([\nabla\_{\hat{\mathbf{x}}} E]\_i)$

Then the time derivative of $E$ is defined as:

$\dot{E} = \left\langle \nabla\_{\hat{\mathbf{x}}} E, \dot{\hat{\mathbf{x}}} \right\rangle = \sum\_i \underbrace{ \langle [\nabla\_{\hat{\mathbf{x}}} E]\_i, \omega_i \rangle}\_{\text{(i)}} - \gamma \underbrace{\langle [\nabla\_{\hat{\mathbf{x}}} E]\_i, \phi\_i([\nabla\_{\hat{\mathbf{x}}} E]\_i) \rangle }\_{\text{(ii)}}$

The $\omega_i$ induces a pure phase rotation; it is orthogonal to the gradient direction, so $\langle[\nabla_{\hat{\mathbf{x}}}E]_i, \omega_i\rangle=0$

And $\phi\_i(\cdot)$ is (i) odd: $\phi_i(-z) = -\phi_i(z)$, and (ii) monotone increasing: $\langle \phi_i(z_1) - \phi_i(z_2), z_1 - z_2 \rangle \geq 0$. With $z = [\nabla_{\hat{\mathbf{x}}} E]_i$, we have $\langle z, \phi_i(z) \rangle \geq 0$, this therefore yields (ii) \leq 0\. Combining (i) and (ii),

$\dot{E}=0-\gamma\sum_i\langle z, \phi_i(z)\rangle\leq0, \quad z=[\nabla_{\hat{\mathbf{x}}} E]_i$

If $W \succeq 0$ and both $w_{ij}$ and $y_i$ are finite, then Eq.(5) is at most quadratic with a finite lower bound:

$E(\hat{\mathbf{x}})\geq-\lambda_{\max}(W) \|\hat{\mathbf{x}}\|^2-\|\mathbf{y}\| \|\hat{\mathbf{x}}\| \geq c_1 \|\hat{\mathbf{x}}\|^2-c_2$, for suitable constants $c_1 > 0$, $c_2 \geq 0$.

Eq.(5) and the bound above provide the classical Lyapunov conditions: $E(\hat{\mathbf{x}}) \geq c_1 \|\hat{\mathbf{x}}\|^2 - c_2, \quad \dot{E}(\hat{\mathbf{x}}) \leq 0$

Hence $E$ is a Lyapunov function for the Eq.(4): it is lower-bounded and monotonically non-increasing along every trajectory, guaranteeing stability.

[Heuristic Analysis: Why BRICK alleviates over-smoothing]

We analyze why the BRICK can theoretically alleviate the oversmoothing by examining its simplified dynamics, steady-state solution and spectral response.

Simplified BRICK Dynamics. To derive an interpretable steady-state solution, we consider a simplified version of BRICK dynamics. We assume a constant natural frequency $\omega_i$ across all nodes and linearize the nonlinear projection $\phi$:

$\frac{d\hat{\mathbf{x}}}{dt} = -\hat{\mathbf{x}}+W\hat{\mathbf{x}}+\mathbf{y}$

This corresponds to a linear consensus-like system, where the $-\hat{\mathbf{x}}$ term mimics a dissipative force, ensuring stability.

Equilibrium Solution. At steady state $( \frac{d\hat{\mathbf{x}}}{dt} = 0)$, we obtain:

$(I - W)\hat{\mathbf{x}}^* = \mathbf{y} \quad \Rightarrow \quad \hat{\mathbf{x}}^* = (I - W)^{-1} \mathbf{y}$

The inverse exists provided the spectral radius of $W$ is strictly smaller than 1 (i.e., $\lambda_k \neq 1, \forall k$), a condition that is typically satisfied for normalized adjacency or attention matrices used in practice.

Spectral Interpretation. Let $W$ be symmetric and decomposed as $W = U \Lambda U^\top$, where $U$ is the orthonormal eigenvector matrix and $\Lambda = \text{diag}(\lambda_1, ..., \lambda_n)$. Projecting into the spectral domain:

$\hat{\mathbf{x}}^* = U(I - \Lambda)^{-1}U^\top \mathbf{y}$

Let $\tilde{\mathbf{y}} = U^\top \mathbf{y}$ and $\tilde{\hat{\mathbf{x}}}^* = U^\top \hat{\mathbf{x}}^*$, then:

$\tilde{\hat{x}}^*_k = \frac{1}{1 - \lambda_k} \tilde{y}_k$

Thus, each spectral mode is scaled by a transfer function: $h(\lambda_k) = \frac{1}{1 - \lambda_k}$

Comparison with Diffusion-based GNNs. In standard diffusion-type GNNs (e.g., GCN), applying $L$ layers is equivalent to using a transfer function $e^{-L\lambda_k}$ in spectral space. This leads to exponential suppression of high-frequency signals (large $\lambda_k$), causing oversmoothing. In contrast, BRICK uses a transfer function that decays much slower: $h(\lambda_k) = \frac{1}{1 - \lambda_k}$, which corresponds to $1/\lambda$-level suppression, allowing high-frequency, discriminative signals to persist.

[1] Rusch, T. K. et al. (2022). Graph-coupled oscillator networks.
[2] Nguyen, T. et al. (2024). From coupled oscillators to graph neural networks: Reducing over-smoothing via a kuramoto model-based approach.
[3] Yeo, B. T. T. et al. (2011). The organization of the human cerebral cortex estimated by intrinsic functional connectivity.
[4] Takeru M. et al. Artificial kuramoto oscillatory neurons.

Thanks again for your feedback, if the reviewer has additional questions, we would be glad to provide further clarifications.

We sincerely thank the reviewer for recognizing the novelty and comprehensive nature of our work. We also greatly appreciate the reviewer’s interest and thoughtful feedback regarding the biological interpretability, computational efficiency and scalability, as well as the optimization dynamics underlying our framework. Below, we address each of the reviewer’s concerns in detail.

1. It would be more impactful to link the synchronization mechanisms to neuroscience domains, and explains why it might be more superior compared to conventional approachs in brain rhythm identification task, and how the discovery in Fig 3 might be relevant to existing studies.

[Link to neuroscience domains]
In principle, the learning behavior of BRICK is shaped by the same governing equations (i.e., Kuramoto model) that describe the neural oscillatory dynamics responsible for generating cognition and behavior. We further introduce structural-functional coupling, implemented through the geometric scattering transform (Line 159), to emulate the cross-frequency coupling observed in cognitive neuroscience. In addition, we incorporate an adaptive control term to enhance task-relevant modulation. Drawing inspiration from neural oscillatory synchronization, our BRICK model not only improves predictive accuracy but also enhances interpretability, thereby offering novel insights into cognitive processes (Section 3.1).

[Advances compared to conventional approaches in brain rhythm identification task]
While traditional models, including CNNs, RNNs, Transformers, and GNNs, have demonstrated strong empirical performance in brain rhythm identification, they largely remain data-driven and agnostic to the underlying neurobiological mechanisms, they often treat the brain as a generic data source without sufficient domain knowledge. In contrast, BRICK explicitly incorporates neural synchronization dynamics, a well-established principle in cognitive neuroscience, into the learning architecture. This alignment with neuroscience allows us to interpret learned patterns (e.g., coupling strengths or synchronization clusters) in biologically meaningful terms (as shown in Figure 3), such as inter-regional communication or pathological desynchronization. Therefore, compared to conventional approaches, BRICK is not just a deep model for feature representation learning, it shows the potential to establish a biologically inspired reasoning system. By embedding core principles of brain dynamics into a learnable, differentiable architecture, our proposed deep model enables achieving more interpretable, temporally-aware, and robust brain state identification.

[Relevance of Figure 3 to existing studies]
Below is a detailed explanation of how the discovered patterns in Figure 3 align with known neurobiological structures and functions:

• Alignment with functional subnetworks (Figure 3a):
  BRICK captures task-specific synchronization across functionally defined brain network in HCP-A. For instance, Visual network regions (red) cluster clearly during the VISMOTOR task in the deeper representation $X^{(L)}$, indicating strong within-network coordination. Similarly, tighter phase clustering emerges in the Sensorimotor and Dorsal Attention networks under tasks such as CART and FACENAME. These results suggest that our model respects functional boundaries in a biologically meaningful manner.
• Intersubject consistency and task-relevant encoding (Figure 3b):
  This figure shows that individuals performing the same task (e.g., “SOCIAL” or “EMOTION” in HCP-YA) form well-separated synchronization in deeper layers (e.g., $X^{(L)}$). This indicates BRICK captures latent dynamics consistent across subjects and datasets, supporting the notion that neural phase alignment underlines cognitive state representation.
• Unsupervised Functional Differentiation (Figure 3c): Here, we assess whether BRICK can produce meaningful cortical parcellations without supervision. Compared to traditional spectral clustering, BRICK yields phase-based partitions that are more spatially localized and better aligned with known functional networks [1]. This supports the claim that BRICK captures biologically coherent dynamics.

By mapping learned feature vectors onto a phase manifold, we can interpret inter-regional phase alignment as a proxy for functional coordination. In this way, BRICK offers a mechanistic perspective grounded in neural oscillation theory, which allows us to move beyond black-box representations toward interpretable latent dynamics.

2. Missing comparisons of the computational cost of the proposed methods compared to existing methods. Would this be more computational expensive compared to conventional graph learning methods, and how it will be scale to larger graph?

[Comparisons of the computational cost]
We have provided a detailed comparison of the inference time per subject for all models across three brain datasets in Appendix A.2 Table 7. In addition, we have now included a summary table that reports the average inference time of each model on all brain datasets. As shown in the results, our proposed model BRICK has slightly higher inference time compared to lightweight baselines such as GCN, GIN, and GraphSAGE. However, it remains more efficient than several other models, while also delivering significantly better performance on brain datasets.

[How to scale to larger graph]
To ensure the scalability of our method to large-scale graphs, we implement all matrix operations using sparse representations. This design allows our model to scale comparably to lightweight models such as GCN, and enables full-batch training on large graphs like ogbn-arxiv without exceeding memory limits.

Moreover, for even larger graphs, e.g., fine-grained applications such as vertex-level cortical modeling in neuroscience, our framework can readily incorporate subgraph sampling strategies [2] to enable smooth training.

3. What might be challenges for the optimization solver for the proposed graph learning procedure, would this guarantee to converge, and how it is sensitive to initialization and hyperparameters?

[Initialization and hyperparameters ablations]
Sensitivity to initialization. To assess the sensitivity of our model to random initialization, we conducted experiments on the Cora, Citeseer, and Pubmed datasets using 5 different random seeds. As reported in Table 2, our model achieves highly consistent results, with the maximum standard deviation across these runs being only 0.36%. This indicates that our method is stable and not significantly affected by initialization variance.

Sensitivity to hyperparameters. We conducted several ablation studies on the most critical hyperparameters and reported them in the appendix.

• In Appendix A.2, we ablate two parameters specific to the BRICK architecture: L (outer iterations for control update) and Q (inner time steps for oscillator dynamics). Increasing both parameters increases inference time, while having only a marginal effect on accuracy. This suggests the model is not overly sensitive to these parameters within a reasonable range.

• In Appendix A.3, we further analyze the sensitivity of geometric scattering transform (GST)-related hyperparameters: level and order. Our results show that increasing either improves performance by capturing more frequency components, but at the cost of slightly increased computation.

[Optimization challenge and converge guarantee]
We appreciate the reviewer’s insightful question regarding the optimization behavior of our proposed method. From a theoretical perspective, the dynamics of our model are governed by a gradient flow is: $\frac{d\hat{x}_i}{dt}=\omega_i+\gamma\cdot \phi_i (y_i+\sum\_{j=1}^{N} w\_{ij}\hat{x}_j)$,

which corresponds to a gradient descent flow on the following energy landscape: $E = - \sum_{i,j} w_{ij} \hat{x}_i^\top \hat{x}_j - \sum_i y_i^\top \hat{x}_i$. This energy function $E$ acts as a natural Lyapunov function, satisfying the condition: $\frac{dE}{dt} \leq 0$, thus ensuring that the system always evolves towards a (local) steady state. Therefore, the entire system remains a valid gradient flow on a constrained manifold, which theoretically guarantees convergence.

In practice, we also observe stable optimization behavior during training. To support this claim, we will include the training loss trajectories of our model on representative datasets in Appendix.

[1] Yeo, B. T. T. et al (2011). The organization of the human cerebral cortex estimated by intrinsic functional connectivity. Journal of Neurophysiology, 106(3), 1125–1165.
[2] Zeng, H. et al (2019). Graphsaint: Graph sampling based inductive learning method. arXiv preprint arXiv:1907.04931.

Thanks again for your feedback, if anything is still unclear or the reviewer would like additional clarification, we would gladly continue the discussion and address any further concerns.

We sincerely thank the reviewer for recognizing the originality of our work and the comprehensiveness of our experimental evaluation. We also appreciate the thoughtful feedback regarding related work and methodology details. Below, we carefully address each concern and provide clarifications or revisions where needed.

1. It is not clear what are the differences between the authors' modified Kuramoto model in BRICK and the ones proposed in [9,10]. What are the advantages of BRICK compared to the work in [44]?

[Differences between our model and [9], [10], [44]]
BRICK and [9]:
While both BRICK and Cabral et al. [9] adopt the Kuramoto model to explore brain dynamics, their roles and applications diverge significantly. In their work, the Kuramoto model is used as a forward simulator to reproduce biologically observed functional connectivity (FC) patterns from empirical structural connectivity (SC), with the goal of explaining emergent brain phenomena. In contrast, our BRICK repurposes the Kuramoto model as a computational module within a machine learning framework. Rather than passively simulating dynamics, BRICK actively learns to perform representation learning, reasoning, and predictive tasks on graphs using dynamic synchronization mechanisms.

BRICK and [10]:
For methodological paradigm, while Capouskova et al. [10] follow a data-driven analytical paradigm, using autoencoders and clustering to uncover latent cognitive states from empirical fMRI data, they encode BOLD phase coherence data using modern machine learning tools without Kuramoto model. Their focus lies in neuroscientific pattern discovery and interpretation. In contrast, BRICK introduces a physics-informed learning framework that embeds the principles of neural synchronization, specifically via the Kuramoto model, directly into the learning process.

The outcomes of these two approaches are also distinct: Capouskova et al. generate empirical insights about cognitive brain states, whereas BRICK yields new machine learning models that not only excel in brain-related tasks but also generalize to broader graph-based AI problems.

BRICK and [44]:
Compared with [44], our BRICK framework offers two key advantages.

• First, neuroscientific grounding: starting from fMRI-derived brain rhythms, we augment the Kuramoto core with a task-driven feedback controller $y_i$ that mimics cognitive control, thereby preserving biological interpretability.

• Second, domain-specific impact: on real HCP fMRI data BRICK achieves the best task decoding and unsupervised parcellation, while on graph benchmarks it delivers competitive accuracy and remains resistant to oversmoothing even at 128 layers.

Taken together, these properties make BRICK a biologically inspired but practically efficient module that bridges brain dynamics and general-purpose AI.

2. Methods. The readability will improve to incorporate more information about the Big-NOS network architecture.

[Details of BIG-NOS]
Conceptually, BIG-NOS can be viewed as an extension of BRICK to graph data. By interpreting graph nodes as analogs of brain regions and the adjacency matrix as a proxy for coupling strength among regions, we seamlessly adapt the physics-informed architecture of BRICK to arbitrary graphs. Both models share the same governing equations rooted in Kuramoto synchronization, but differ in their data domains and learning objectives.

In BRICK, the inputs consist of BOLD signals over brain regions along with structural or functional connectivity; in BIG-NOS, the inputs are general node features and graph topology. Likewise, the outputs diverge in form: BRICK predicts cognitive or disease-related brain states (analogous to graph-level classification), while BIG-NOS predicts node/graph labels in standard graph learning settings.

We are committed to include the formal network detail of BIG-NOS into Section 3.2 as part of the overall architecture description.

3. Can you clarify why equation 4 can be interpreted as gradient flow?

[Gradient flow interpretation of Eq.(4)]

Gradient flow refers to a class of dynamical systems that evolve along the steepest descent direction of an energy functional $E$. The evolution equation takes the form $\dot{\hat{x}}_i = -\frac{\partial E}{\partial \hat{x}_i}$

which guarantees that the energy $\dot{E}$ decreases over time, i.e., $\dot{E} = \frac{dE}{dt} \leq 0$

leading the system toward a local minimum of $E$. The energy function defined in Eq.(5) $E = -\sum_{i,j} w_{ij} \hat{x}_i^\top \hat{x}_j - \sum_i y_i^\top \hat{x}_i$ encodes two key objectives:

• The first term promotes pairwise alignment (synchronization) among oscillators $\hat{x}\_i$ and $\hat{x}\_j$ with strong coupling $w\_{ij}$.

• The second term aligns each oscillator with its task-specific control pattern $y_i$ (controlled synchronization).

Taking the gradient of $E$ with respect to $\hat{x}_i$ yields $\frac{\partial E}{\partial \hat{x}\_i} = -\sum_j w\_{ij} \hat{x}\_j - y\_i$,

so that the negative gradient becomes $-\frac{\partial E}{\partial \hat{x}\_i} = \sum_j w\_{ij} \hat{x}\_j + y_i$

Substituting this into Eq.(4) we obtain: $\dot{\hat{x}}_i = \omega_i + \gamma \phi_i \left( -\frac{\partial E}{\partial \hat{x}_i} \right)$

Thus, Eq.(4) is considered a gradient flow because its dynamics are equivalent to evolving along the negative gradient direction of an energy functional $E$. Here, the operator $\phi$ does not destroy the gradient flow structure, but instead imposes structural constraints on the descent direction (e.g., ensuring motion remains on the unit sphere).

4. What kind on nerural nertwork architecture did the author use to parameterize the natural frequency Ω (lines 201-202)?

[Parameterization of natural frequency Ω]
The parametrization of natural frequency Ω does not rely on a deep neural network. Instead, Ω is parameterized using a set of learnable vectors whose norms determine the rotation speed (i.e., natural frequency) of oscillators. In the forward pass, these frequency values are used to rotate the 2D oscillator states (like turning a point in a circle), simulating how each node evolves over time due to its intrinsic dynamics.

[9] Joana Cabral, Etienne Hugues, Olaf Sporns, and Gustavo Deco. Role of local network oscillations in resting-state functional connectivity. Neuroimage, 57(1):130–139, 2011.373
[10] Katerina Capouskova, Morten L Kringelbach, and Gustavo Deco. Modes of cognition: Evidence from metastable brain dynamics. NeuroImage, 260:119489, 2022.375
[44] Takeru Miyato, Sindy Löwe, Andreas Geiger, and Max Welling. Artificial kuramoto oscillatory neurons. arXiv preprint arXiv:2410.13821, 2024

Should there remain any ambiguities or if the reviewer has further questions, we would be more than happy to engage in continued discussion and resolve any remaining issues.

We sincerely thank the reviewer for acknowledging the novelty and the value of our work. Below, we provide detailed responses to all concerns.

[W1 poor exp setting]，[W5 unsupported claims]

[Exp setting and model selection]

We acknowledge that some of the reported accuracy for baseline models are lower than those in prior studies that specifically optimized each model/dataset. However, our intention was to ensure a fair and controlled comparison across all models under uniform exp settings, particularly important consideration in neuroscience, where the focus lies in uncovering mechanistic and functional insights.

For all models, we fix major hyperparams as hidden dim 256; #layers 4 (2 for GCN/GAT); #epochs 1500; lr 5e-4 – 1e-3; weight decay 5e-4. We will report all the details in final draft.

While these settings may not be optimal for every model, e.g., GCNII perform best with 32–64 layers on Cora,which is not applicable for GCN or GIN due to over-smoothing or instability. Our choice thus ensures comparability under identical architectural constraints. We would also like to emphasize that the accuracy of GCN on Cora in our setting (81.66%) is highly consistent with prior reports (e.g., 81.5% in GRAND), which validates our implementation.

In addition, since our main focus is on neuroscience-related applications aligns with the Neuroscience and cognitive science track, we prioritize baseline models that have been validated in computational neuroscience, such as GCN, GAT, GIN etc. [1]. Nevertheless, we have extended our comparisons by including GTN (2019), GRAND (2021), GraphCON (2022) and KuramotoGNN (KGNN) (2024). The results show that our models achieve significant improvements in brain datasets and retain strong performance in general graph benchmarks.

Brain states identification:

Node classification:

Graph classification:

[W2 lack of method justification]

We now clarify the motivation, justification, and necessity of our method. Our primary focus is brain rhythm identification from neural oscillations, which are fundamental in coordinating communication and cognitive functions [2]. Moreover, Kuramoto model has long been used to study coupled neural oscillators [3]. Building on these foundation, BRICK explicitly models oscillatory dynamics of brain function by Kuramoto dynamics, which not only achieves strong performance in brain state decoding (Tab.2), but also yields interpretable, biologically meaningful synchronization patterns (Fig.3). We believe this justifies both the necessity and novelty of our approach in the perspective of computational neuroscience.

Encouraged by the promising results on brain data, we sought to genelize neural oscillation principles to general graph domains in AI/ML, which share the same notion of coupled dynamical systems as in brain networks. We hypothesize that learning distinctive graph representations is analogous to coordinating functional roles (graph-level) across brain regions (node-level).

We also provide the following methodology justifications. If there are specific aspects of the theoretical formulation that the reviewer would like further clarification on, we would be more than happy to provide detailed explanations as needed.

Proof of Lyapunov conditions.
Let $\hat{\mathbf{x}} = [\hat{x}\_1, ..., \hat{x}\_N] \in \mathbb{R}^{N \times d}$ and recall dynamical system Eq.(4) and energy $E$ Eq.(5). Because $W = [w\_{ij}]$ is symmetric (undirected),

$\nabla_{\hat{\mathbf{x}}} E = - (W \hat{\mathbf{x}} + \mathbf{y}), \quad \mathbf{y} = [y_1 \ldots y_N]$

Define $u_i := y_i + \sum_j w\_{ij} \hat{x}\_j$. Then Eq.(4) may be rewritten as:

$\dot{\hat{x}}_i = \omega_i + \gamma \phi\_i(u_i) = \omega\_i - \gamma \phi\_i(-u\_i) = \omega\_i - \gamma \phi\_i([\nabla\_{\hat{\mathbf{x}}} E]\_i)$

Then the time derivative of $E$ is defined as:

$\dot{E} = \left\langle \nabla\_{\hat{\mathbf{x}}} E, \dot{\hat{\mathbf{x}}} \right\rangle = \sum\_i \underbrace{ \langle [\nabla\_{\hat{\mathbf{x}}} E]\_i, \omega_i \rangle}\_{\text{(i)}} - \gamma \underbrace{\langle [\nabla\_{\hat{\mathbf{x}}} E]\_i, \phi\_i([\nabla\_{\hat{\mathbf{x}}} E]\_i) \rangle }\_{\text{(ii)}}$

The $\omega_i$ induces a pure phase rotation; it is orthogonal to the gradient direction, so $\langle[\nabla_{\hat{\mathbf{x}}}E]_i, \omega_i\rangle=0$

And $\phi\_i(\cdot)$ is (i) odd: $\phi_i(-z) = -\phi_i(z)$, and (ii) monotone increasing: $\langle \phi_i(z_1)-\phi_i(z_2), z_1-z_2\rangle\geq 0$. With $z = [\nabla_{\hat{\mathbf{x}}} E]_i$, we have $\langle z, \phi_i(z) \rangle \geq 0$, this therefore yields (ii) \leq 0\. Combining (i) and (ii),

$\dot{E}=0-\gamma\sum_i\langle z, \phi_i(z)\rangle\leq0, \quad z=[\nabla_{\hat{\mathbf{x}}} E]_i$

If $W \succeq 0$ and both $w_{ij}$ and $y_i$ are finite, then Eq.(5) is at most quadratic with a finite lower bound:

$E(\hat{\mathbf{x}})\geq-\lambda_{\max}(W) \|\hat{\mathbf{x}}\|^2-\|\mathbf{y}\| \|\hat{\mathbf{x}}\| \geq c_1 \|\hat{\mathbf{x}}\|^2-c_2$, for suitable constants $c_1 > 0$, $c_2 \geq 0$.

Eq.(5) and the bound above provide the classical Lyapunov conditions: $E(\hat{\mathbf{x}}) \geq c_1 \|\hat{\mathbf{x}}\|^2 - c_2, \quad \dot{E}(\hat{\mathbf{x}}) \leq 0$

Hence $E$ is a Lyapunov function for the Eq.(4): it is lower-bounded and monotonically non-increasing along every trajectory, guaranteeing stability.

[W3 missing details]

[Details of BRICK and BIG-NOS]

The governing equation presents in Tab.1 for BIG-NOS is shared with BRICK. This alignment is intentional and central to our design, as it enables a seamless adaptation of BRICK from human brain to general graph learning tasks. We'll give more details in final version.

[Interpretation of Fig.3]

We now elaborate how the discovered patterns in Fig. 3 support neurobiological plausibility:

• Alignment with functional subnetworks (Fig.3a):
  BRICK captures task-specific synchronization across functional brain network in HCP-A (e.g., Visual network regions (red) cluster clearly during the VISMOTOR task in the deeper feature $X^{(L)}$, indicating strong within-network coordination). It suggest that BRICK respects functional boundaries.
• Inter-subject consistency and task-relevant encoding (Fig.3b):
  Subjects performing the same task (e.g., “SOCIAL” in HCP-YA) form well-separated synchronization in deeper $X^{(L)}$. This indicates BRICK captures consistent latent dynamics across subjects, supporting that neural phase alignment underlines cognitive state representation.
• Unsupervised functional differentiation (Fig.3c): BRICK can produce meaningful cortical parcellations without supervision. Compared to spectral clustering, BRICK yields phase-based partitions that are more spatially aligned with known functional networks [4].

To provide stronger quantitative evidence, we have already included numeric results in Tab.2 and Appendix A.4. However, we will also include the following numeric figures to support Fig.3:

• Within-network vs. between-network phase variance for Fig.3a to quantify the degree of synchronization.
• Purity and NMI to measure how well the unsupervised parcels in Fig.3c align with the existing functional network.

[W4 missing prior works]

We have carefully reviewed both papers and will include them in the paper with a discussion of differences.

1. Compare with GraphCON: GraphCON is based on second-order ODEs representing damped nonlinear oscillators. Its formulation is inspired by mechanical oscillators to simulate oscillatory dynamics. BRICK is explicitly informed by neural oscillation, which is first-order and phase-based, aiming to model functional coordination in brain dynamics and translate those principles into learnable graph architectures.

2. Compare with KuramotoGNN: Their work demonstrates the use of synchronization dynamics for improving message passing in GNN, thus mitigating oversmoothing. However, BRICK incorporates vector-valued oscillators and adaptive control mechanisms. we ground our model in neuroscientific applications and demonstrate its performance and interpretability on neural data, beyond generic graph benchmarks.

[1] Dan, T. et al. (2024). Exploring the enigma of neural dynamics through a scattering-transform mixer landscape for riemannian manifold.
[2] Fries, P. (2005). A mechanism for cognitive dynamics: Neuronal communication through neuronal coherence.
[3] Cabral, J.et al. (2011). Role of local network oscillations in resting-state functional connectivity.
[4] Yeo, B.T. et al. (2011). The organization of the human cerebral cortex estimated by intrinsic functional connectivity.

We are glad to provide further clarification or discuss in any format that would be helpful.

Dear authors,

Thank you for your rebuttal. I have carefully read your rebuttal and would like to raise few more points and questions.

[On W1 poor experimental setting]. First, I am still not convinced about the fairness and rigor of experimental setting.

• I agree that some baseline performance (e.g., GCN) seem reasonable.
• on fairness: However, the authors used unreasonably detrimental hyperparameters for GCNII (specifically, small number of layers), and thus, I do not consider it a fair comparison. If the authors want to focus on the specific hyperpameter space, I strongly encourage the authors to remove the baselines that were not intended to be used with the hyperparameters. Likewise, I think the updated experiment with the added baselines is misleading as well, as both GraphCon and Grand have been reported to perform best at layers deeper than the ones the authors used.
• on the focus of the paper: While the authors claimed that their focus is on “neuroscience-related applications aligns with the Neuroscience and cognitive science track”, they reported more results about graph learning (2 Tables and 1 Figure; 12 benchmark datasets) than those of neuroscience applications (1 Table and 1 Figure; 3 benchmark datasets). Moreover, in lines 314-315, the authors claim that their method “achieves SOTA performance on both heterophilic and homophilic graphs”. With such a heavy focus on graph learning, I must insist that the authors should not sidestep from the issue I raised (W1 poor experimental setting) by claiming they are not the focus. Besides, the reported performances of the proposed method is far away from SOTA performance. I encourage the authors to avoid using misleading claims.

[On W2 lack of method justification].

• I understand that the authors designed their method by leveraging Kuramoto model and neural oscillation, which understandably may improve “brain learning”.
• However, the authors’ motivation goes beyond brain learning. Specifically, they aimed to “rethink graph learning through the lens of coupled neural oscillators” (line 58) and design “a new graph learning mechanism via neural oscillatory synchronization” (line 90).
• Why is leveraging neural oscillation and kuramoto model important for “graph learning”? As I stated earlier, the only motivation that the authors provided is speculative (”each node acts as a coupled oscillator, evolving through interactions governed by the graph topology”; line 78-79). In summary, the link between proposed method’s design principles and graph learning is unclear and unjustified.
• Thus, contrary to the authors’ claim, I am not convinced this method is important for graph learning.

In the future version of the paper, I encourage the authors to either (1) substantially reduce their focus on graph learning or (2) implement more rigorous and standard experimental settings accepted in graph learning domain, with a clear demonstration of the superiorities introduced by the proposed method over existing graph learning methods. In the current version, with such a heavy focus on graph learning and its poor justification/experiments, I cannot recommend acceptance.

Sincerely,

We thank the reviewer for the constructive feedback and for considering a score adjustment based on the conceptual value of our work. Below, we provide a more focused response to your concerns.

1. However, I generally do not consider 'promising results' to be good enough for a top-tier publication, especially with only thin experiments on brain data.

we appreciate the opportunity to clarify the breadth and depth of our experimental validation on brain data.

In the current submission, we included three brain task datasets (HCP-A, HCP-YA, and HCP-WM) where we quantitatively demonstrate our model’s performance (Table 2), along with visualizations of synchronized neural dynamics in feature space (Figure 3a and 3b). We further explored the model’s ability to uncover biologically meaningful structures without supervision for functional parcellation (Figure 3c). Together, these results aim to show both quantitative and qualitative merits of our method in brain applications.

Beyond the main manuscript, we have in fact conducted more extensive evaluations on neurological disease datasets, which were omitted due to space constraints, and we just focused on task-evoked fMRI datasets. We now summarize the results for three disease-related datasets (ADNI, PPMI, NIFD) here to provide a broader perspective:

• Alzheimer’s Disease Neuroimaging Initiative (ADNI): includes 135 resting-state fMRI samples from subjects diagnosed with Alzheimer’s disease (AD) or cognitively normal (CN) controls.

• Parkinson’s Progression Markers Initiative (PPMI): includes 173 samples spanning Parkinson’s disease (PD), SWEDD (scans without evidence of dopaminergic deficit), prodromal, and healthy controls.

• Neuroimaging Initiative for Frontotemporal Lobar Degeneration (NIFD): includes 1010 samples across a spectrum of frontotemporal dementia subtypes, including logopenic variant of primary progressive aphasia (LSD), behavioral variant (BV), progressive non-fluent aphasia (PNFA), semantic variant (SV), and cognitively normal (CON).

On all three datasets, our model maintains strong classification performance. In particular, the results on PPMI and NIFD are statistically significant (*, $p < 0.01$, paired t-test), confirming the robustness of the observed improvements. This highlights the model’s potential in handling complex, heterogeneous brain disorders.

We believe this is because BRICK captures abnormal coordination patterns between brain regions, which is often a hallmark of neurodegeneration. The latent phase representation enables the model to detect subtle disruptions in functional connectivity that might be missed by conventional GNNs.

We hope this addresses the reviewer’s concern by demonstrating that BRICK is consistently effective across a wide spectrum of brain datasets.