Uncover Governing Law of Pathology Propagation Mechanism Through A Mean-Field Game

We sincerely thank the reviewer for their thoughtful and constructive comments. We truly appreciate the recognition of our work’s strengths, including the innovative integration of mean-field game theory, physics-informed modeling, and symbolic regression for modeling tau–Aβ interactions in Alzheimer's disease. We are encouraged by the reviewer’s positive evaluation of the model’s interpretability, methodological clarity, and superior predictive performance compared to a broad spectrum of existing approaches. Below, we address each of the reviewer’s concerns point by point.

Q1: Fitting an independent reaction law at each of the 100K surface vertices greatly increases computation and memory costs.

A1: Thank you for raising this important concern. We agree that fitting an independent reaction law at ~100,000 surface vertices is computationally intensive. To address this, we use a GPU-efficient symbolic regression model that enforces shared functional forms with sparse parameterization. Dropout-based sparsification further enhances convergence and reduces overfitting.

Training is conducted on an 8×NVIDIA H100 NVL GPU system (94GB per GPU), enabling parallel computation over large-scale data. While training takes $\sim$10h, inference is highly efficient ($\sim$0.2s). Despite the computational demands, this vertex-level modeling is essential to capture the regional heterogeneity of tau–Aβ interactions, which would otherwise be lost in coarser, atlas-based approaches. We will clarify these details in the final version.

Q2: The surface - based mesh graph constructed from cortical triangular tessellation is well - suited for simulating local trans - neuronal spread, but it may have limitations in capturing long - range propagation.

A2: Thank you for the insightful suggestion. We acknowledge that the current surface-based mesh primarily captures local, cortical trans-neuronal spread. To address long-range propagation mechanisms,we are currently developing a multi-layer diffusion framework that integrates both local (surface-based) and global (connectome-based) long-range pathways. This hybrid approach allows us to simulate tau spread across multiple spatial scales, reflecting both short-range anatomical adjacency and long-range anatomical or functional connectivity. This layered formulation is aligned with recent neurobiological models and will be a key component of our ongoing and future work.

Q3: How reliable is the symbolic regression component in inferring biologically meaningful tau–Aβ reaction laws?

A3: Thank you for your helpful comment. We agree that validation against known biology is essential. To address this, we now provide an explicit interpretation of a representative learned symbolic equation for tau update in the medial parahippocampal region. The expression includes five interpretable terms:

• $u_{44}$ (S_calcarine): The negative coefficient from the calcarine cortex (a primary visual region affected in late Braak stages V–VI) may reflect a dampening influence on tau accumulation in early-stage regions—potentially capturing long-range regulatory effects in network dynamics [Braak & Braak, 1991].
• $v_{81}$ (S_cingul-Mid-Ant): The positive linear Aβ contribution from the anterior cingulate cortex—a region known to accumulate Aβ early (Thal Phase I–II)—suggests a modest but biologically plausible driving force for tau spread [Thal et al., 2002; Palmqvist et al., 2017].
• $\sin(2.2 \cdot v_9)$ (G_cingul-Post-dorsal): A nonlinear oscillatory term, potentially reflecting a thresholded or biphasic Aβ influence, aligns with evidence that Aβ’s impact on tau propagation may saturate or even reverse at high concentrations [Busche & Hyman, 2020].
• Constant term (+1.9): Captures the intrinsic baseline accumulation of tau in the medial temporal lobe, consistent with spontaneous age-related tauopathy in entorhinal and parahippocampal areas [Braak & Del Tredici, 2015].
• $\frac{1.9}{\exp(3.5 \cdot v_{16}) + 1}$ (G_front_sup–parahippocampal): A sigmoid-shaped Aβ term, indicative of a saturating dose-response, reflects known dynamics where tau is more sensitive to Aβ at subthreshold levels but less responsive once amyloid burden becomes extensive [Jack et al., 2013].

These regions collectively form the entorhinal–cingulate–frontal circuit, known as the epicenter of early tau pathology (Braak I–III). Importantly, these symbolic terms were not cherry-picked; they consistently emerged as top contributors across folds (see Fig. 6), providing converging evidence that our model captures biologically grounded and interpretable tau–Aβ interaction rules.

We will add this analysis to the final version.

Refs:

[Braak & Braak, 1991] Neuropathological staging of Alzheimer-related changes, Acta Neuropathologica

[Thal et al., 2002] Phases of Aβ-deposition in the human brain and its relevance for the development of AD, Neurology

[Palmqvist et al., 2017] Earliest accumulation of β-amyloid occurs within the default-mode network and concurrently affects brain Synergy between amyloid-β and tau in Alzheimer's disease connectivity, Nature Communications

[Busche & Hyman, 2020] Synergy between amyloid-β and tau in Alzheimer's disease, Nature Neuroscience

[Jack et al., 2013] Tracking pathophysiological processes in Alzheimer's disease: an updated hypothetical model of dynamic biomarkers, Lancet Neurology

[Braak & Del Tredici, 2015] The preclinical phase of the pathological process underlying sporadic Alzheimer's disease

Q4: Is the model’s performance limited by the quality or resolution of the neuroimaging datasets?

A4: This is a good question. The performance of our model does depend on the resolution and quality of the neuroimaging data, as tau and amyloid SUVR values are derived from PET images with limited spatial resolution and potential noise. To mitigate these issues, we model tau propagation on the cortical surface mesh, which preserves anatomical geometry and supports fine-grained modeling of local diffusion and interregional interactions (Fig. 2a). Furthermore, by integrating symbolic regression and connectivity-aware diffusion (via GNNs), our model can partially compensate for noise by learning consistent propagation dynamics across folds and cohorts.

Additionally, our GAN-based optimal transport formulation operates on tau density distributions, rather than raw voxel-wise intensities, making it more robust to imaging variability. As shown in Table 1, MFG4AD achieves strong predictive performance across both ADNI and OASIS, despite differences in scanner protocols and cohort demographics. Future work may explore incorporating higher-resolution imaging modalities or harmonization techniques to further enhance reliability. Thanks.

We sincerely thank the reviewer for their thoughtful and encouraging feedback. We are especially grateful for your recognition of the novelty of our formulation, the use of realistic cortical geometry, and the biological interpretability introduced by the symbolic regression module. Your positive evaluation motivates us to further refine and extend this line of research. Below, we address each of your concerns in detail.

Q1: The model assumes that amyloid influences tau, but doesn’t explore the reverse, or model amyloid propagation at all.

A1: Thank you for raising this important point. Our model focuses on Aβ → tau interactions because substantial neuropathological and longitudinal imaging studies suggest that amyloid-β (Aβ) accumulation precedes and facilitates tau pathology, rather than the reverse. Aβ and tau exhibit distinct propagation mechanisms: Aβ spreads in a relatively diffuse and spatially global manner—primarily governed by extracellular deposition—while tau propagates transneuronally in a more localized and network-constrained fashion (Braak et al., 1991; Thal et al., 2002; Vogel et al., 2020).

Moreover, Aβ accumulation tends to plateau early in the disease course, while tau continues to spread and correlate more strongly with neurodegeneration and cognitive decline (Jack et al., 2013; Hanseeuw et al., 2019). Given these asymmetries in dynamics and temporal ordering, we adopt the prevailing neuroscience perspective by modeling Aβ as an upstream driver of tau, not the other way around.

Future extensions could incorporate Aβ propagation using different mechanisms (e.g., global diffusion with clearance dynamics), but doing so would require introducing distinct priors or differential equations, which would make the model significantly more complex. Our current scope prioritizes biological plausibility and model interpretability for tau dynamics, where Aβ acts as an external input shaping regional vulnerability. We will add this discussion to the final version.

Q2: The learned symbolic terms are interesting, but could benefit from more validation or comparison to known biology.

A2: Thank you for your very helpful comment. We agree that validation against known biology is crucial. To address this, we now provide an explicit interpretation of a representative learned symbolic equation for tau update at the medial parahippocampal region:

• $u_{44}$ (S_calcarine): The negative coefficient from the calcarine cortex (a primary visual region affected in late Braak stages V–VI) may reflect a dampening influence on tau accumulation in early-stage regions—potentially capturing long-range regulatory effects in network dynamics [Braak & Braak, 1991].
• $v_{81}$ (S_cingul-Mid-Ant): The positive linear Aβ contribution from the anterior cingulate cortex—a region known to accumulate Aβ early (Thal Phase I–II)—suggests a modest but biologically plausible driving force for tau spread [Thal et al., 2002; Palmqvist et al., 2017].
• $\sin(2.2 \cdot v_9)$ (G_cingul-Post-dorsal): A nonlinear oscillatory term, potentially reflecting a thresholded or biphasic Aβ influence, aligns with evidence that Aβ’s impact on tau propagation may saturate or even reverse at high concentrations [Busche & Hyman, 2020].
• Constant term (+1.9): Captures the intrinsic baseline accumulation of tau in the medial temporal lobe, consistent with spontaneous age-related tauopathy in entorhinal and parahippocampal areas [Braak & Del Tredici, 2015].
• $\frac{1.9}{\exp(3.5 \cdot v_{16}) + 1}$ (G_front_sup–parahippocampal): A sigmoid-shaped Aβ term, indicative of a saturating dose-response, reflects known dynamics where tau is more sensitive to Aβ at subthreshold levels but less responsive once amyloid burden becomes extensive [Jack et al., 2013].

These regions collectively form the entorhinal–cingulate–frontal circuit, known as the epicenter of early tau pathology (Braak I–III). Importantly, these symbolic terms were not cherry-picked; they consistently emerged as top contributors across folds (see Fig. 6), providing converging evidence that our model captures biologically grounded and interpretable tau–Aβ interaction rules.

We will add this analysis into the final version.

Q3: The paper doesn’t address uncertainty in predictions, which could be important in clinical settings.

A3: This is a great comment. We agree with the reviewer that incorporating uncertainty estimates could be valuable for clinical applications. While our current model provides point estimates, we plan to extend it with Bayesian inference or dropout-based uncertainty quantification in future work. Thanks.

Q4: A difference map would also be informative.

A4: Thank you for the suggestion. We have added the difference maps for both the ADNI and OASIS datasets (see Fig. R2). These maps highlight regional deviations between predicted and observed tau levels, offering additional insight into model performance. We have prepared the corresponding population-level brain mapping (Rebuttal_Fig.R2.pdf) to support verification. However, due to the rebuttal policy restricting supplementary uploads during the review process, we are unable to provide this figure at this stage. Nevertheless, we will incorporate these results into the final version of the manuscript, as we agree that they provide important support for the paper. We sincerely appreciate your suggestion.

Q5: In general, along with the weaknesses mentioned above, the lack of temporal analyses would add much value to this work. Also, I am wondering if there could be more subject-level analyses to appreciate the quality of the individualized predictions.

A5: Thank you for this valuable suggestion. In our response, we have added subject-level (16 subjects across different disease stages on ADNI dataset) scatter plots comparing predicted versus observed tau values on ADNI datasets (see Fig. R3, will be public after acceptance), providing a clearer view of individualized prediction quality.

Additionally, we agree that incorporating temporal analyses would further enrich the work. While this current study focuses on baseline-to-follow-up prediction, we are actively extending our framework to model full longitudinal trajectories, which we plan to include in future work.

Q6: Mean Filed Game -> Mean Field Game

A6: Amended. Thanks.

Refs:

[1] Braak, H., & Braak, E. (1991). Neuropathological stageing of Alzheimer-related changes. Acta neuropathologica.

[2] Thal, D. R., et al. (2002). Phases of Aβ-deposition in the human brain and its relevance for the development of AD. Neurology

[3] Palmqvist, S., et al. (2017). Earliest accumulation of β-amyloid occurs within the default-mode network and concurrently affects brain connectivity. Nature communications

[4] Busche, M. A., & Hyman, B. T. (2020). Synergy between amyloid-β and tau in Alzheimer’s disease. Nature neuroscience [5] Braak, H., & Del Tredici, K. (2015). The preclinical phase of the pathological process underlying sporadic Alzheimer’s disease. Brain

[6] Jack, C. R., et al. (2013). Tracking pathophysiological processes in Alzheimer's disease: an updated hypothetical model of dynamic biomarkers. The lancet neurology

Dear Reviewer 5WWy,

As the deadline for author-reviewer discussion is approaching, could you please check the authors' rebuttal and post your response?

Thank you!

Best,

AC

We thank the reviewer for recognizing our work. MFG4AD combines physics-based modeling and deep learning in a novel way to track tau spread directly on the brain surface. It not only predicts future tau accurately but also reveals interpretable Aβ–tau interactions, showing strong results on ADNI and OASIS datasets. Below, we address each of your concerns in detail.

Q1: limited direct qualitative or quantitative evidence to substantiate how interpretable the learned laws

A1: Thank you for highlighting the need for more explicit evidence supporting the interpretability of the learned Aβ–tau interaction law. In response, we have revised Section 3.2 (line 298-320) to provide a detailed examination of the symbolic expressions $R_j(u,v)$ generated by our model.

To demonstrate mechanistic interpretability, we present a representative symbolic reaction equation at a vertex in the medial parahippocampal cortex—a region known to initiate tau pathology (Braak I–II) [Braak & Braak, 1991]. The symbolic expression incorporates terms from spatially and functionally linked Aβ sources, including the anterior/posterior cingulate cortex and frontal–parahippocampal junction, consistent with known structural pathways and disease staging trajectories. Each term in the symbolic expression is biologically meaningful:

• Positive coefficients (e.g., $+0.02*v_{81}$​) denote facilitatory effects from local or upstream Aβ burden. For instance, anterior cingulate cortex ($v_{81}$) has been identified as an early Aβ hub in Thal Phase 2 and is tightly linked to the medial temporal lobe via the cingulum bundle [Thal et al., 2002,Neurology; Jacobs et al., 2018,Neuron].

• Nonlinear components (e.g., $0.2\sin(2.2v_9) + 1.9 + 1.9\exp(3.5v_{16}) + 1$) capture thresholded or saturating influences of Aβ, reflecting empirical findings that Aβ’s effect on tau becomes pronounced only beyond certain levels [Hanseeuw et al., 2019, JAMA Neurology]. These forms also simulate sigmoid-like molecular switch behavior, analogous to pathological tipping points in disease cascades.

• Negative weights (e.g., $-0.09u_{44}$​) arising from regions such as the calcarine (visual cortex) suggest inhibitory or damping effects, consistent with their involvement in later Braak stages (V–VI) and possible negative feedback or non-contributory roles in early tau spread [Busche & Hyman, 2020,Nature Neuroscience].

Importantly, the spatial configuration of contributing terms reflects a biologically coherent propagation route—the entorhinal–hippocampal–cingulate loop, widely regarded as the epicenter of early Alzheimer’s pathology [Braak & Braak, 1991; Jacobs et al., 2018].

Together, these findings demonstrate that the learned symbolic laws are not only mathematically predictive, but also biologically interpretable, aligning with canonical disease progression stages and known anatomical connectivity.

Q2: insufficient technical details regarding its symbolic regression module

A2 Thank you for this important point. We have provided more detailed implementation of the symbolic regression module. Specifically:

• Architecture: The symbolic module (SymbolicNetL00) is composed of multiple SymbolicLayerL0 layers (as detailed in the anonymous GitHub). Each layer transforms inputs using a collection of interpretable mathematical functions (e.g., identity, sin, exp, product) and learns to select sparse combinations through a differentiable binary mask.

• Search space of functions: The model uses a library of symbolic primitives including: Unary functions: identity, $\sin(x)$, $\cos(x)$, $\exp(x)$, etc. Binary functions: product $x \cdot y$, quotient $x/y$, etc.

• Optimization & regularization: We employ the hard concrete distribution as in [Louizos et al., 2017] to sample differentiable binary masks $z \in [0,1]$ applied element-wise to the weight matrix. This enables $L0$ sparsity regularization, encouraging simpler and more interpretable symbolic formulas.

• Selection criteria: After training, symbolic terms with high dropout probabilities (i.e., low expected activation z) are pruned to remove weak or unstable components. The remaining active terms are ranked based on the magnitude of their learned weights, and their corresponding expressions are reconstructed using a predefined set of interpretable mathematical functions (e.g., linear, sine, exponential), ensuring biological and mathematical interpretability.

• Robustness & Reproducibility: We evaluated the learned symbolic expressions using 5-fold cross-validation. Key contributing regions—such as the entorhinal cortex—consistently ranked among the top terms across folds. Fig. 6 summarizes the selection frequency of each Aβ region, highlighting the stability and biological plausibility of the learned interaction patterns. We will make it clear in the final version, thanks.

Q3: the sensitivity of the model's performance to variations in these hyperparameters

A3: Thank you for your comment, we have included the sensitivity analysis of hyperparameter choices (line 203-208 in the manuscript).

• A smaller λ fails to adequately enforce reconstruction fidelity, leading to suboptimal predictions. The best performance is observed when $λ \in [10, 15]$, indicating this range strikes a desirable balance between adversarial realism and numerical accuracy.
• Gradient penalty (gp) shows mild sensitivity across 5–15, suggesting the 1-Lipschitz constraint is sufficiently enforced within this range without overly penalizing the critic’s learning.
• Critic steps ($n_C$) shows no monotonic trend; moderate values (e.g., 5) suffice to balance critic convergence and training efficiency.
• Learning rates ($\eta_C = 1\text{e-5}, \eta_G = 1\text{e-4}$) yield the best result, as they support stable adversarial updates. Higher rates destabilize optimization, while lower rates may slow convergence.

Training time is 74.02s/it, inference time is 0.23 s, on H100 NVL with 8 GPUs (94 G/per).

Q4: the biological assumptions about tau propagation, and empirical evidence or analysis show how the derived reaction laws

A4: Thank you for the question. Our model treats each cortical vertex as an agent whose tau level evolves under local interactions, governed by a surface-based brain graph. This graph is constructed directly from the native cortical surface mesh, preserving geodesic continuity and the anatomical folding pattern of the cortex (see Appendix A.1, Fig. 7). Unlike SC/FC-based connectomes, our geometry-aware network ensures that tau propagation follows biologically realistic, anatomically constrained trajectories.

The weak coupling assumption in MFG aligns well with the biological reality: tau spread in one region is influenced not by a single source but by the average tau level across its neighboring regions. This naturally captures the notion of local accumulation and distributed influence along the cortical mantle. We adapt the MFG framework by integrating spatial connectivity into the interaction terms, such that the propagation dynamics respect surface-based neighborhood structures.

The learned reaction laws provide not only enhanced predictive capability but also interpretable symbolic expressions, which allows us to quantify how local Aβ burden and spatially adjacent influences drive tau accumulation. These symbolic laws highlight early-stage epicenters (e.g., entorhinal and medial prefrontal cortex) and map out downstream propagation patterns consistent with Braak I–III staging. By coupling biophysically grounded diffusion, geometric network constraints, and interpretable equations, our model offers novel, data-driven insights into the spatial logic of tau spread in Alzheimer's disease.

Q5: lacks a comprehensive discussion

A5: Thank you for this important question. Beyond predictive accuracy, our model generates interpretable diffusion fields that offer biological insights into tau propagation dynamics. In Fig. 5, we visualize the spatial propagation of tau across clinical stages and datasets (ADNI and OASIS). The inferred flow fields $\nu$ reveal coherent directional patterns that mirror canonical staging: for instance, early spread originates in the medial temporal lobe (e.g., entorhinal cortex), progressing toward isocortical regions such as parietal and frontal cortices—consistent with Braak stages I–VI.

Importantly, these flow fields were derived purely from data via our model, without imposing priors about direction or staging. The emergence of biologically plausible propagation axes, repeatedly observed across datasets and stages, demonstrates that our framework captures not just the distribution but also the mechanism of tau spread. This provides novel empirical support for network-based propagation hypotheses in AD.

Dear Reviewer mzBH,

We would like to sincerely thank you for the time, effort, and expertise you devoted to reviewing our submission. Your thoughtful feedback, particularly regarding the hyperparameter setting analysis and interpretability of the learned laws, has helped us significantly improve the quality of our work. We are truly grateful for the opportunity to improve our work through this revision process.

We hope that our rebuttal has adequately addressed your concerns. If so, we would be most grateful if you would consider updating your score accordingly.

Of course, if you have any remaining concerns or further suggestions, please do not hesitate to let us know, we would be more than happy to address them to the best of our ability.

With sincere appreciation,

Authors

Thank you for your thoughtful summary and evaluation of our work. We greatly appreciate your recognition of our framework’s strengths—including its principled integration of biophysical modeling and deep learning, the novel MFG-based formulation, and its potential to yield both accurate predictions and interpretable mechanistic insights. Below, we address each of your concerns in detail.

Q1: Subgroups' performance

A1: We appreciate the reviewer’s attention to subgroup-level performance. In the original manuscript, we provided visual comparisons of model performance and tau propagation across diagnostic subgroups (CN, EMCI, LMCI, and AD), as shown in Fig. 5 (main text) and Fig. 8 (Appendix). Table 1 summarizes the performance on all samples. Thank you for reminding us and this is a very helpful suggestion. We have now included the quantitative performance metrics (e.g., MAE, RMSE) for each subgroup. These results complement the original table and provide a more rigorous and transparent evaluation of model performance across disease stages. We will incorporate all the results in the final version.

We also note that the CN group in our analysis includes subjects with Subjective Memory Complaint (SMC) in ADNI dataset. While clinically categorized as cognitively normal, SMC individuals often present subtle biomarker changes or early cognitive concerns. Their inclusion likely increases the heterogeneity within the CN group and may contribute to the slightly higher prediction error observed in this subgroup. We have clarified this point in the revised manuscript.

Q2: Relevant literature about Abeta-Tau interactions and, quantitative cross‑dataset evidence

A2: We appreciate this observation and have revised the manuscript to provide both (i) supporting literature and (ii) quantitative cross‑dataset evidence.

• We now cite several large‑scale imaging and modelling studies demonstrating that local Aβ burden amplifies connectivity‑constrained tau propagation: [1] Lee et al. “Regional A-tau interactions promote onset and accelertion of Alzheimer’s disease tau spreading”, Neuron, 2022. [2] Roemer-Casiano et al. “Amyloid-associated hyperconnectivity drives tau spread across connected brain regions in Alzheimer’s disease”, Science Translatioal Medicine, 2025
• To demonstrate cross‑dataset consistency, we computed vertex‑wise Aβ–tau interaction frequency maps for ADNI (64,262 left‑ and 54,128 right‑hemisphere vertices) and OASIS (34,917 LH / 26,898 RH) and then summed these frequencies within the 160 Destrieux regions of interest [Destrieux, C., et al. Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature. Neuroimage]. This procedure identified 67 high‑frequency ROIs in ADNI and 37 in OASIS; 22 regions (24663 vertices)—bilateral entorhinal, parahippocampal, fusiform, inferior/middle temporal, posterior cingulate/precuneus, caudate, and amygdala—were common in both cohorts, encompassing more than 20 000 vertices (Dice = 0.42). Crucially, the Aβ‑related frequencies in these 22 overlapping parcels remained strongly concordant across datasets (Spearman ρ = 0.46). A 10,000‑iteration permutation test that randomly reassigned OASIS ROI labels produced an empirical p‑value of 0.009, confirming that such concordance is exceedingly unlikely to occur by chance and reinforcing the cross‑dataset reproducibility of the spatial pattern of Aβ–tau synergy observed in ADNI and OASIS.

We will incorporate all the results in the final version.

Q3: The correlation between the truth and the predicted one

A3: Thank you for your comment. We have added a dedicated figure (Fig. R1) that visualises predicted vs. observed tau‑SUVR means for every diagnostic subgroup in each cohort. Each dot is an individual subject, the numerical results are summarised here:

Predicted vs. Observed Tau-SUVR Fitting Metrics

These metrics demonstrate an agreement between predicted and true values across all subgroups. We have prepared the corresponding fitting curves (Rebuttal_Fig.R1.pdf) to support verification. However, due to the rebuttal policy restricting supplementary uploads during the review process, we are unable to provide this figure at this stage. Nevertheless, we will incorporate these results into the final version of the manuscript, as we agree that they provide important support for the paper. We sincerely appreciate your suggestion.

Thank you for your time and consideration.