Re-Think and Re-Design Graph Neural Networks in Spaces of Continuous
  Graph Diffusion Functionals
Graph neural networks (GNNs) are widely used in domains like social networks
and biological systems. However, the locality assumption of GNNs, which limits
information exchange to neighboring nodes, hampers their ability to capture
long-range dependencies and global patterns in graphs. To address this, we
propose a new inductive bias based on variational analysis, drawing inspiration
from the Brachistochrone problem. Our framework establishes a mapping between
discrete GNN models and continuous diffusion functionals. This enables the
design of application-specific objective functions in the continuous domain and
the construction of discrete deep models with mathematical guarantees. To
tackle over-smoothing in GNNs, we analyze the existing layer-by-layer graph
embedding models and identify that they are equivalent to l2-norm integral
functionals of graph gradients, which cause over-smoothing. Similar to
edge-preserving filters in image denoising, we introduce total variation (TV)
to align the graph diffusion pattern with global community topologies.
Additionally, we devise a selective mechanism to address the trade-off between
model depth and over-smoothing, which can be easily integrated into existing
GNNs. Furthermore, we propose a novel generative adversarial network (GAN) that
predicts spreading flows in graphs through a neural transport equation. To
mitigate vanishing flows, we customize the objective function to minimize
transportation within each community while maximizing inter-community flows.
Our GNN models achieve state-of-the-art (SOTA) performance on popular graph
learning benchmarks such as Cora, Citeseer, and Pubmed.
Tingting Dan, Jiaqi Ding, Ziquan Wei, Shahar Z Kovalsky, Minjeong Kim, Won Hwa Kim, Guorong Wu
# Re-Think and Re-Design Graph Neural Networks in Spaces of Continuous Graph Diffusion Functionals

###### Abstract

Graphs are ubiquitous in various domains, such as social networks and biological systems. Despite the great successes of graph neural networks (GNNs) in modeling and analyzing complex graph data, the inductive bias of locality assumption, which involves exchanging information only within neighboring connected nodes, restricts GNNs in capturing long-range dependencies and global patterns in graphs. Inspired by the classic _Brachistochrone_ problem, we seek how to devise a new inductive bias for cutting-edge graph application and present a general framework through the lens of variational analysis. The backbone of our framework is a two-way mapping between the discrete GNN model and continuous diffusion functional, which allows us to design application-specific objective function in the continuous domain and engineer discrete deep model with mathematical guarantees. _First_, we address over-smoothing in current GNNs. Specifically, our inference reveals that the existing layer-by-layer models of graph embedding learning are equivalent to a \(\ell_{2}\)-norm integral functional of graph gradients, which is the underlying cause of the over-smoothing problem. Similar to edge-preserving filters in image denoising, we introduce the total variation (TV) to promote alignment of the graph diffusion pattern with the global information present in community topologies. On top of this, we devise a new selective mechanism for inductive bias that can be easily integrated into existing GNNs and effectively address the trade-off between model depth and over-smoothing. _Second_, we devise a novel generative adversarial network (GAN) to predict the spreading flows in the graph through a neural transport equation. To avoid the potential issue of vanishing flows, we tailor the objective function to minimize the transportation within each community while maximizing the inter-community flows. Our new GNN models achieve state-of-the-art (SOTA) performance on graph learning benchmarks such as _Cora_, _Citeseer_, and _Pubmed_.

## 1 Introduction

Graph is a fundamental data structure that arises in various domains, including social network analysis [48], natural language processing [42], computer vision [34], recommender systems [43], and knowledge graphs [25] among others. Tremendous efforts have been made to operate machine learning on graph data (called graph neural networks, or GNNs) at the node [27], link [47], and graph level [31]. The common inductive bias used in GNNs is the homophily assumption that nodes that are connected in a graph are more likely to have similar features or labels. In this context, most GNN models deploy a collection of fully-connected layers to progressively learn graph embeddings by
aggregating the nodal feature representations from its topologically-connected neighbors throughout the graph [21].

Under the hood of GNNs, the graph representation learning process is achieved by various learnable operations, such as message passing [5] or graph convolution [27]. Due to the nature of exchanging information in a local graph neighborhood, however, it is challenging to capture global graph representations, which go beyond node-to-node relationship, by leveraging the deep architecture in GNNs while being free of overly smoothing the feature representations for the closely-connected nodes. Fig. 1 demonstrates the root cause of over-smoothing issue in current GNNs, where node color denotes the group label (no color means unlabeled) and edge thickness indicates connection strength. It is clear that nodes #1 and #2 are located at the boundary of two communities. The inductive bias of GNNs (i.e., locality assumption) enforces the node embedding vectors on node #1 and #2 becoming similar due to being strongly connected (highlighted in red), even though the insight of global topology suggests that their node embeddings should be distinct. As additional layers are added to GNNs, the node embeddings become capable of capturing global feature representations that underlie the entire graph topology. However, this comes at the cost of over-smoothing node embeddings across graph nodes due to (1) an increased number of node-to-node information exchanges, and (2) a greater degree of common topology within larger graph neighborhoods. In this regard, current GNNs only deploy a few layers (typically two or three) [30], which might be insufficient to characterize the complex feature representations on the graph.

It is evident that mitigating the over-smoothing problem in GNNs will enable training deeper models. From a network architecture perspective, skip connections [18; 45], residual connections [29; 20], and graph attention mechanisms [39; 37] have been proposed to alleviate the information loss in GNNs, by either preserving the local feature representation or making information exchange adaptive to the importance of nodes in the graph. Although these techniques are effective to patch the over-smoothing issue in some applications, the lack of an in-depth understanding of the root cause of the problem poses the challenge of finding a generalized solution that can be scaled up to current graph learning applications.

Inspired by the success of neural ordinary differential equations in computer vision [10], research focus has recently shifted to link the discrete model in GNNs with partial differential equation (PDE) based numerical recipes [44; 33; 6; 15]. For example, Graph Neural Diffusion (GRAND) formulates GNNs as a continuous diffusion process [6]. In their framework, the layer structure of GNNs corresponds to a specific discretization choice of temporal operators. Since PDE-based model does not revolutionize the underlying inductive bias in current GNNs, it is still unable to prevent the excessive information change between adjacent nodes as in nodes #1 and #2 in Fig. 1. In this regard, using more advanced PDE solvers only can provide marginal improvements in terms of numerical stability over the corresponding discrete GNN models, while the additional computational cost, even in the feed-forward scenario, could limit the practical applicability of PDE-based methods for large-scale graph learning tasks.

In this regard, pioneering work on continuous approaches has prompted to re-think GNN as a graph diffusion process governed by the Euler-Lagrange (E-L) equation of the heat kernel. This formulation is reminiscent of the _Brachistochrone_ problem 1, which emerged over 400 years ago and established the mathematical framework of classical mechanics. The powerful calculus of variations allows us to generate solutions for various mechanics questions (e.g., the slope that yields the fastest ball sliding down the curve is given by a cycloid) through the lens of E-L equation, as shown in Fig. 2 (top).

Figure 1: Demonstration of the root cause of over-smoothing in GNNs. Nodes #1 and #2 are located along the boundary of two communities. The locality assumption in GNNs steers the learning of the graph representations by constraining the information exchange via node-to-node connections. However, such link-wise inductive bias opts to neutralize the contrast of node embeddings between nodes #1 and #2, which might undermine the node classification accuracy. Our research framework yields the solution for the over-smoothing issue by enabling heat-kernel diffusion within each community while penalizing the excessive community-to-community information exchanges (highlighted in red).
In a similar vein, the question that arises in the context of community detection is: What graph diffusion pattern is best suited for preserving community organizations? The question for graph classification would be: What graph diffusion pattern works best for capturing the system-level characteristics of graph topology? Following the spirit of _Brachistochrone_ problem, we present a general research framework to customize application-specific GNNs in a continuous space of graph diffusion functionals. As shown in Fig. 2 (bottom), we have established a fundamental structure for our framework that involves a two-way mapping between a discrete GNN model and a continuous graph diffusion functional. This allows us to develop application-specific objective functions (with an explainable regularization term) in the continuous domain and construct a discrete deep model with mathematical guarantee. We demonstrate two novel GNN models, one for addressing over-smoothing and one for predicting the flows from longitudinal nodal features, both achieving state-of-the-art performance (_Cora_: 85.6%, _Citeseer_: 73.9%, _Pubmed_: 80.10%, even in 128 network layers).

We have made four major contributions. (1) We establish a connection between the discrete model of GNNs and the continuous functional of inductive bias in graph learning by using the E-L equation as a stepping stone to bridge the discrete and continuous domains. (2) We introduce a general framework to re-think and re-design new GNNs that is less "black-box". (3) We devise a novel selective mechanism upon inductive bias to address the over-smoothing issue in current GNNs and achieve state-of-the-art performance on graph learning benchmarks. (4) We construct a novel GNN in the form of a generative adversarial network (GAN) to predict the flow dynamics in the graph by a neural transport equation.

## 2 Methods

In the following, we first elucidate the relationship between GNN, PDE, and calculus of variations (COV), which sets the stage for the _GNN-PDE-COV_ framework for new GNN models in Section 2.2.

Re-think GNNs: Connecting dots across graph neural networks, graph diffusion process, Euler-Lagrange equation, and Lagrangian mechanics

**Graph diffusion process.** Given graph data \(\mathcal{G}=(V,W)\) with \(N\) nodes \(V=\left\{v_{i}|i=1,\ldots,N\right\}\), the adjacency matrix \(W=\left[w_{ij}\right]_{i,j=1}^{N}\in\mathbb{R}^{N\times N}\) describes connectivity strength between any two nodes. For each node \(v_{i}\), we have a graph embedding vector \(x_{i}\in\mathcal{R}^{m}\). In the context of graph topology, the graph gradient \(\left(\nabla_{\mathcal{G}}x\right)_{ij}=w_{ij}\left(x_{i}-x_{j}\right)\) indicates the feature difference between \(v_{i}\) and \(v_{j}\) weighted by the connectivity strength \(w_{ij}\), where \(\nabla_{\mathcal{G}}\) is a \(\mathbb{R}^{N}\rightarrow\mathbb{R}^{N\times N}\) operator. Thus, the graph diffusion process can be formulated as \(\frac{\partial x(t)}{\partial t}=div\left(\nabla_{\mathcal{G}}x(t)\right)\), where the evolution of embedding vectors \(x=\left[x_{i}\right]_{i=1}^{N}\) is steered by the graph divergence operator.

**Connecting GNN to graph diffusion.** In the regime of GNN, the regularization in the loss function often measures the smoothness of embeddings \(x\) over the graph by \(x^{T}\Delta x\), where \(\Delta=div(\nabla_{\mathcal{G}})\) is the graph Laplacian operator [27]. To that end, the graph smoothness penalty encourages two connected nodes to have similar embeddings by information exchange in each GNN layer. Specifically, the new graph embedding \(x^{l}\) in the \(l^{th}\) layer is essentially a linear combination of the graph embedding \(x^{l-1}\) in the previous layer, i.e., \(x^{l}=A_{W,\Theta}x^{l-1}\), where the matrix \(A\) depends on graph adjacency matrix

Figure 2: Top: The _Brachistochrone_ problem played a pivotal role in the development of classical mechanics and the powerful mathematical tool known as the calculus of variations. Bottom: We introduce a general framework to answer _”Brachistochrone”_ problems regarding diffusion patterns on graphs that allows us to re-think and re-design application-specific deep model of GNNs with enhanced mathematical interpretability.
\(W\) and trainable GNN parameter \(\Theta\). After rewriting \(x^{l}=Ax^{l-1}\) into \(x^{l}-x^{l-1}=(A-I)x^{l-1}\), updating graph embeddings \(x\) in GNN falls into a discrete graph diffusion process, where the time parameter \(t\) acts as a continuous analog of the layers in the spirit of Neural ODEs [10]. It has been shown in [6] that running the graph diffusion process for multiple iterations is equivalent to applying a GNN layer multiple times.

**GNN is a discrete model of Lagrangian mechanics via E-L equation.** The diffusion process \(\frac{\partial x(t)}{\partial t}=div\left(\nabla_{\mathcal{G}}x(t)\right)\) has been heavily studied in image processing in decades ago, which is the E-L equation of the functional \(\min\limits_{x}\int_{\Omega}|\nabla x|^{2}dx\). By replacing the 1D gradient operator defined in the Euclidean space \(\Omega\) with the graph gradient \(\left(\nabla_{\mathcal{G}}x\right)_{ij}\), it is straightforward to find that the equation governing the graph diffusion process \(\frac{\partial x(t)}{\partial t}=div\left(\nabla_{\mathcal{G}}x(t)\right)\) is the E-L equation of the functional \(\min\limits_{x}\int_{\mathcal{G}}|\nabla_{\mathcal{G}}x|^{2}dx\) over the graph topology. Since the heat kernel diffusion is essentially the mathematical description of the inductive bias in current GNNs, we have established a mapping between the mechanics of GNN models and the functional of graph diffusion patterns in a continuous domain.

**Tracing the smoking gun of over-smoothing in GNNs.** In Fig. 1, we observed that the inductive bias of link-wise propagation is the major reason for excessive information exchange, which attributes to the over-smoothing problem in GNNs. An intuitive approach is to align the diffusion process with high-level properties associated with graph topology, such as network communities. After connecting the GNN inductive bias to the functional of graph diffusion process, we postulate that the root cause of over-smoothing is the isotropic regularization mechanism encoded by the \(\ell_{2}\)-norm. More importantly, connecting GNN to the calculus of variations offers a more principled way to design new deep models with mathematics guarantees and model mechanistic explainability.

### Re-design GNNs: Revolutionize inductive bias, derive new E-L equation, and construct deeper GNN

The general roadmap for re-designing GNNs typically involves three major steps: (1) formulating inductive bias into the functional of graph diffusion patterns; (2) deriving the corresponding E-L equation; and then (3) developing a new deep model of GNN based on the finite difference solution of E-L equation. Since the graph diffusion functional is application-specific, we demonstrate the construction of new GNN models in the following two graph learning applications.

#### 2.2.1 Develop _VERY_ deep GNNs with a selective mechanism for link-adaptive inductive bias

**Problem formulation.** Taking the feature learning component (learnable parameters \(\Theta\)) out of GNNs, the graph embeddings \(x^{L}\) (output of an \(L\)-layer GNN) can be regarded as the output of an iterative smoothing process (\(L\) times) underlying the graph topology \(\mathcal{G}\), constrained by the data fidelity \(\left\|x^{L}-x^{0}\right\|_{2}^{2}\) (w.r.t. the initial graph embeddings \(x^{0}\)) and graph smoothness term \(\int_{\mathcal{G}}|\nabla_{\mathcal{G}}x|^{2}dx\). Inspired by the great success of total variation (TV) for preserving edges in image denoising [35], reconstruction [41] and restoration [8], we proposed to address the over-smoothing issue in current GNN by replacing the quadratic Laplacian regularizer with TV on graph gradients, i.e., \(\mathcal{J}_{TV}(x)=\int|\nabla_{\mathcal{G}}x|dx\). Thus, the TV-based objective function for graph diffusion becomes: \(\min\limits_{x}\bigl{(}\left\|x-x^{0}\right\|_{2}^{2}+\mathcal{J}_{TV}(x))\).

However, the \(\ell_{1}\)-norm function, denoted by \(|\cdot|\) in the definition of the total variation functional \(\mathcal{J}_{TV}\), is not differentiable at zero. Following the dual-optimization schema [4; 7], we introduce the latent auxiliary matrix \(z\in\mathbb{R}^{N\times N}\) and reformulate the TV-based functional as \(\mathcal{J}_{TV}(x,z)=\max\limits_{z}\min\limits_{x}\int(z\otimes\nabla_{ \mathcal{G}}x)dx\), subject to \(|z|\leq\mathbf{1}^{N\times N}\), where \(\otimes\) is Hadamard operation between two matrices. Furthermore, we use an engineering trick of element-wise operation \(z_{ij}(\nabla_{\mathcal{G}}x)_{ij}\) to keep the degree always non-negative (same as we take the absolute value), which makes the problem solvable. In the end, we reformulate the minimization of \(\mathcal{J}_{TV}(x)\) into a dual _min-max_ functional as \(\mathcal{J}_{TV}(x,z)\), where we maximize \(z\) (\(z\to\mathbf{1}^{N\times N}\)) such that \(\mathcal{J}_{TV}(x,z)\) is close enough to \(\mathcal{J}_{TV}(x)\). Therefore, the new objective function is reformulated as:

\[\mathcal{J}(x,z)=\max\limits_{z}\min\limits_{x}\left\|x-x^{0}\right\|_{2}^{2}+ \lambda\int(z\nabla_{\mathcal{G}}x)dx, \tag{1}\]
which \(\lambda\) is a scalar balancing the data fidelity term and regularization term. Essentially, Eq. 1 is the dual formulation with _min-max_ property for the TV distillation problem [50].

**Constructing E-L equations.** To solve Eq. 1, we present the following two-step alternating optimization schema. _First_, the inner minimization problem (solving for \(x_{i}\)) in Eq. 1 can be solved by letting \(\frac{\partial}{\partial x_{i}}\mathcal{J}(x_{i},z_{i})=0\):

\[\frac{\partial}{\partial x_{i}}\mathcal{J}(x_{i},z_{i})=2(x_{i}-x_{i}^{0})+ \lambda z_{i}\nabla_{\mathcal{G}}x_{i}=0\ \ \Rightarrow\ \ \hat{x}_{i}=x_{i}^{0}-\tfrac{\lambda}{2}z_{i}\nabla_{\mathcal{G}}x_{i} \tag{2}\]

Replacing \(\left(\nabla_{\mathcal{G}}x\right)_{ij}\) with \(w_{ij}\left(x_{i}-x_{j}\right)\), the intuition of Eq. 2 is that each element in \(\hat{x}_{i}\) is essentially the combination between the corresponding initial value in \(x_{i}^{0}\) and the overall graph gradients \(z_{i}\nabla_{\mathcal{G}}x_{i}=\sum_{j\in\mathcal{N}_{i}}w_{ij}(x_{i}-x_{j}) z_{i}\) within its graph neighborhood \(\mathcal{N}_{i}\). In this regard, Eq. 2 characterizes the dynamic information exchange on the graph, which is not only steered by graph topology but also moderated by the attenuation factor \(z_{i}\) at each node.

_Second_, by substituting Eq. 2 back into Eq. 1, the objective function of \(z_{i}\) becomes \(\mathcal{J}(z_{i})=\max\limits_{|x_{i}|\leq 1}\left\|\tfrac{\lambda}{2}z_{i} \nabla_{\mathcal{G}}x_{i}\right\|_{2}^{2}+\lambda z_{i}\nabla_{\mathcal{G}}(x _{i}^{0}-\tfrac{\lambda}{2}z_{i}\nabla_{\mathcal{G}}x_{i})\). With simplification (in Eq. 8 to Eq. 10 of Supplementary), the optimization of each \(z_{i}\) is achieved by \(\operatorname*{arg\,min}\limits_{|z_{i}|\leq 1}z_{i}\nabla_{\mathcal{G}}x_{i}z _{i}\nabla_{\mathcal{G}}x_{i}-\tfrac{4}{\lambda}z_{i}\nabla_{\mathcal{G}}x_{i} ^{0}\). Specifically, we employ the majorization-minimization (MM) method [16] to optimize \(z_{i}\) by solving this constrained minimization problem (the detailed derivation process is given in Eq. 8 to Eq. 19 of Section 5.1 of Supplementary), where \(z_{i}\) can be iteratively refined by:

\[z_{i}^{l}=clip\underbrace{(z_{i}^{l-1}+\frac{2}{\beta\lambda}\nabla_{ \mathcal{G}}x_{i},1)}_{b}=\left\{\begin{array}{cc}b&|b|\leq 1\\ 1&b>1\\ -1&b<-1\end{array}\right. \tag{3}\]

\(\beta\) is a hyper-parameter that is required to be no less than the largest eigenvalue of \((\nabla_{\mathcal{G}}x_{i})(\nabla_{\mathcal{G}}x_{i})^{\intercal}\).

**Develop new GNN network architecture with a selective inductive bias.** The building block in vanilla GNN [27] is a FC (fully-connected) layer where the input is the embedding vectors after isotropic graph diffusion (in \(\ell_{2}\)-norm). Since the estimation of graph embeddings \(x\) in Eq. 2 depends on the latest estimation of \(z^{(l)}\), such recursive _min-max_ solution for Eq. 1 allows us to devise a new network architecture that disentangles the building block in vanilla GNN into the feature representation learning and graph diffusion underling TV. As shown in Fig. 3, we first deploy a FC layer to update the graph embeddings \(x^{(l)}\). After that, we concatenate a diffusion-clip (DC) layer for selective graph diffusion, which sequentially applies (1) node-adaptive graph diffusion (blue arrow in Fig. 3) on \(x^{(l)}\) by Eq. 22, and (2) clip operation (purple arrow in Fig. 3) to each \(x_{i}^{(l)}\) by Eq. 3.

Footnote 2: Since the optimization schema has been switched to the layer-by-layer manner, the initialization \(x_{0}\) becomes \(x^{(l-1)}\) from the previous layer.

**Remarks.** Eq. 3 indicates that larger connective degree results in larger value of \(z\). Thus, the DC layer shifts the diffusion patterns by penalizing the inter-community information exchange (due to strong connections) while remaining the heat-kernel diffusion within the community. The preference of such link-adaptive diffusion can be adjusted by the hyper-parameter \(\lambda\)3 in Eq. 1. Recall our intuitive solution for over-smoothing problem in Fig. 1, the DC layer offers the exact global insight of graph topology to keep the node embeddings distinct between nodes #1 and #2. We demonstrate the effect of DC layer on the real-world graph data in Fig. 8 of Supplementary document.

Figure 3: Our new deep model integrates a novel diffusion-clip (DC) layer (for selective graph diffusion) after the conventional fully-connected (FC)layer (for graph representation learning).
#### 2.2.2 Predict flow dynamics through graph neural transport equation

**Problem formulation.** We live in a world of complex systems, where everything is intricately connected in multiple ways. A holistic insight of how the system's components interact with each other and how changes in one part of the system can affect the behavior of the whole sheds new light on the dynamic behaviors of these complex systems over time. However, oftentimes it is an ill-posed problem. Taking the toy system in Fig. 4(a) as an example, while it is simple to calculate the future focal patterns based on the focal patterns at the current time point and the node-to-node flow information, determining flow dynamics based on longitudinal nodal observations is computationally hard since the solution is not unique.

The naive solution to predict the spreading flow is to (1) train a GNN to learn the intrinsic node embeddings and (2) predict the flow based on the difference of learned embeddings. However, this two-step approach might suffer from vanishing flow due to over-smoothing in GNNs. Following the spirit of _Brachistochrone_ problem, we ask the question "What flow field \(f(t)=[f_{ij}(t)]_{i,j=1}^{N}\) underlines the system mechanics to the extent that it is able to predict the behaviors in the future?"

In this paper, we focus on the conservative system of energy transportation [2]. The system mechanics is formulated as:

\[\frac{dx}{dt}+div(q)=0 \tag{4}\]

where \(q=[q_{ij}]_{i,j=1}^{N}\) is the flux field which propagates the potential energy \(u(t)=[u_{i}(t)]_{i=1}^{N}\) (conserved quantity) over time. Similar to a gravity field driving water flow, the intuition of Eq. 4 is that the change of energy density \(u\) (we assume there is a non-linear mapping \(\phi\) from external force \(x\) to \(u\), i.e., \(u_{i}=\phi(x_{i})\)) leads to energy transport throughout the entire graph. As flux is closely related to the difference of energy \(\nabla_{\mathcal{G}}u\) underlying the graph topology, we assume the energy flux \(q\) is regulated by the potential energy field \(\nabla_{\mathcal{G}}u\), i.e., \(q=\alpha\otimes\nabla_{\mathcal{G}}u\), where \(\alpha=[\alpha_{ij}]_{i,j=1}^{N}\) is a learnable matrix characterizing the link-wise contribution of each energy potential \(\nabla_{\mathcal{G}}u_{ij}\) to the potential energy flux \(q_{ij}\).

By plugging \(q=\alpha\otimes\nabla_{\mathcal{G}}u\) into Eq. 4, the energy transport process can be reformulated as:

\[\frac{\partial u}{\partial t}=-\phi^{-1}(\alpha\otimes div(\nabla_{\mathcal{G }}u))=-\phi^{-1}(\alpha\otimes\Delta u), \tag{5}\]

where \(\Delta=div(\nabla_{\mathcal{G}})\) is the graph Laplacian operator. Since the PDE in Eq. 5 is equivalent to the E-L equation of the quadratic functional \(\mathcal{J}(u)=\min\limits_{u}\int_{\mathcal{G}}\alpha\otimes|\nabla_{ \mathcal{G}}u|^{2}du\) (after taking \(\phi\) away), a major issue is the over-smoothness in \(u\) that might result in vanishing flows. In this context, we propose to replace the \(\ell_{2}\)-norm integral functional \(\mathcal{J}(u)\) with TV-based counterpart \(\mathcal{J}_{TV}(u)=\min\limits_{u}\int_{\mathcal{G}}\alpha\otimes|\nabla_{ \mathcal{G}}u|du\).

**Renovate new E-L equation - graph neural transport equations.** Following the _min-max_ optimization schema in Eq. 1-3, we introduce an auxiliary matrix \(f\) to lift the undifferential-able barrier. After that, the minimization of \(\mathcal{J}_{TV}(u)\) boils down into a dual _min-max_ functional \(\mathcal{J}_{TV}(u,f)=\min\limits_{u}\max\limits_{f}\int_{\mathcal{G}}\alpha \otimes f(\nabla_{\mathcal{G}}u)du\). Since \(u(t)\) is a time series, it is difficult to derive the deterministic solutions (as Eq. 2-3) by letting \(\frac{\partial}{\partial u}\mathcal{J}_{TV}=0\) and \(\frac{\partial}{\partial f}\mathcal{J}_{TV}=0\). Instead, we use \(G\ddot{a}\acute{e}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute {a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a}\acute{a} \acute{a}\acute
**Remarks.** The solution to Eq. 6 is known as continuous max-flow and constitutes a continuous version of a graph-cut [1]. Since \(\alpha\) is a latent variable and potential energy \(u\) is given, the maximization of \(f\) opts towards maximizing the spreading flow through the lens of \(\alpha\). In this regard, the mechanistic role of auxiliary matrix \(f\) is essentially the latent (maximized) spreading flows that satisfy \(u(t+1)_{i}=u(t)_{i}+\sum_{j=1}^{N}f_{ij}(t)\). The potential energy \(\hat{u}\) can be solved via a wave equation (\(u_{tt}=div(f_{t})=\alpha^{2}\otimes\Delta u\)), where the system dynamics is predominated by the adjusted Lagrangian mechanics \(\alpha^{2}\otimes\Delta u\). By determining \(\alpha\) at a granularity of graph links, we devise a novel GAN model to predict the spreading flows \(f\) which not only offers explainability underlying the _min-max_ optimization mechanism in Eq. 6 but also sets the stage to understand system dynamics through machine learning.

**Develop a GAN model of flow prediction with TV-based Lagrangian Mechanics.** The overall network architecture is shown in Fig. 4 (b), which consists of a generator (red solid box) and a discriminator module (blue solid box). Specifically, the generator (\(G\)) consists of (1) a GCN component [15] to optimize \(\hat{u}\) through the wave equation \(u_{tt}=\alpha^{2}\otimes\Delta u\) and (2) a FC layer to characterize the non-linear mapping function \(\hat{x}(t+1)=\phi^{-1}(\hat{u}(t))\). In contrast, the discriminator (\(D\)) is designed to (1) synthesize \(\alpha\) and (2) construct the future \(\tilde{u}_{t+1}\) based on the current \(u_{t}\) and current estimation of spreading flow \(f=\alpha\otimes\nabla_{G}u\) (orange dash box). To make the network architecture consistent between generator and discriminator modules, we include another GCN to map the synthesized \(\tilde{u}(t+1)\) to the external force \(\tilde{x}(t+1)\). Note, since the working mechanism of this adversarial model underlines the min-max optimization in the energy transport equation, the nature of predicted spreading flows is carved by the characteristics of _max-flow_.

The driving force of our network is to minimize (1) the MSE (mean square error) between the output of the generator \(\hat{x}_{t+1}\) and the observed regional features, (2) the distance between the synthesized regional features \(\tilde{x}_{t+1}\) (from the discriminator) and the output of generator \(\hat{x}_{t+1}\) (from the generator). In the spirit of probabilistic GAN [49], we use one loss function \(\mathcal{L}_{D}\) to train the discriminator (\(D\)) and another one \(\mathcal{L}_{G}\) to train the generator (\(G\)):

\[\left\{\begin{array}{c}\mathcal{L}_{D}=D\left(x_{t+1}\right)+\left[\xi-D \left(G\left(x_{t}\right)\right)\right]^{+}\\ \mathcal{L}_{G}=D\left(G\left(x_{t}\right)\right)\end{array}\right. \tag{7}\]

where \(\xi\) denotes the positive margin and the operator \([\cdot]^{+}=\max(0,\cdot)\). Minimizing \(\mathcal{L}_{G}\) is similar to maximizing the second term of \(\mathcal{L}_{D}\) except for the non-zero gradient when \(D(G(x_{t}))\geq\xi\).

## 3 Experiments

In this section, we evaluate the performance of the proposed _GNN-PDE-COV_ framework with comparison to six graph learning benchmark methods on a wide variety of open graph datasets [36], as well as a proof-of-concept application of uncovering the propagation pathway of pathological events in Alzheimer's disease (AD) from the longitudinal neuroimages.

### Datasets and experimental setup

**Dataset and benchmark methods.** We evaluate the new GNN models derived from our proposed GNN framework in two different applications. _First_, we use three standard citation networks, namely _Cora_, _Citeseer_, and _Pubmed_[36] for node classification (the detailed data statistic is shown in Table 3 of Supplementary). We adopt the public fixed split [46] to separate these datasets into training, validation, and test sets. We follow the experimental setup of [9] for a fair comparison with six benchmark GNN models (vanilla GCN [27], GAT [39], GCNII [9], ResGCN [29], DenseGCN [29], GRAND [6]). Since our DC-layer can be seamlessly integrated into existing GNNs as a plug-in. The corresponding new GNN models (with DC-layer) are denoted GCN+, GAT+, GCNII+, ResGCN+, DenseGCN+, and GRAND+, respectively.

_Second_, we apply the GAN model in Section 2.2.2 to predict the concentration level of AD-related pathological burdens and their spreading pathways from longitudinal neuroimages. Currently, there is no _in-vivo_ imaging technique that can directly measure the flow of information across brain regions. Here, our computational approach holds great clinical value to understand the pathophysiological mechanism involved in disease progression [26]. Specifically, we parcellate each brain into 148 cortical surface regions and 12 sub-cortical regions using Destrieux atlas [13]. The wiring topology of these 160 brain regions is measured by diffusion-weighted imaging [3] and tractography techniques
[19]. The regional concentration levels AD pathology including amyloid, tau, and fluorodeoxyglucose (FDG) and cortical thickness (CoTh) are measured from PET (positron emission tomography) and MRI (magnetic resonance imaging) scans [23]. We use a total of \(M=1,291\) subjects from ADNI [32], each having longitudinal imaging data (2-5 time points). The details of image statistics and pre-processing are shown in Sec. 6.1.2. Since we apply the flow prediction model to each modality separately, we differentiate them with _X-FlowNet_ (\(X\) stands for amyloid, tau, FGD, and CoTh).

**Experimental setup**. In the node classification task, we verify the effectiveness and generality of DC-layer in various number of layers (\(L=2,4,8,16,32,64,128\)). All baselines use their own default parameter settings. Evaluation metrics include accuracy, precision and F1-score. To validate the performance of _X-FlowNet_, we examine (1) prediction accuracy (MAE) of follow-up concentration level, (2) prediction of the risk of developing AD using the baseline scan, and (3) understand the propagation mechanism in AD by revealing the node-to-node spreading flows of neuropathologies.

The main results of graph node classification and flow prediction are demonstrated in Section 3.2 and 3.3, respectively. Other supporting results such as ablation study and hyper-parameter setting are shown in Section 6 of the Supplementary document.

\begin{table}
\begin{tabular}{c l l l l l l l l} \hline \hline
**Dataset** & **Model** & \(L=2\) & \(L=4\) & \(L=8\) & \(L=16\) & \(L=32\) & \(L=64\) & \(L=128\) \\ \hline \multirow{4}{*}{**Cora**} & GCN & 81.30 & 79.90 & 75.70 & 25.20 & 20.00 & 31.80 & 20.80 \\  & **GCN+** & **82.761\({}^{*}_{*,007}\)** & **82.726\({}^{*}_{*,007}\)** & **82.306\({}^{*}_{*,007}\)** & **70.604\({}^{*}_{*,04}\)** & **67.804\({}^{*}_{*,08}\)** & **66.605\({}^{*}_{*,04}\)** & **59.903\({}^{*}_{*,017}\)** \\ \cline{2-10}  & GAT & 82.40 & 80.30 & 57.90 & 31.90 & - & - & - \\  & **GAT+** & **82.600\({}_{*,020}\)** & **80.500\({}_{*,020}\)** & **69.701\({}^{*}_{*,01}\)** & **66.003\({}^{*}_{*,01}\)** & **63.60\({}^{*}_{*,01}\)** & **54.60\({}^{*}_{*,04}\)** & **45.70\({}^{*}_{*,07}\)** \\ \cline{2-10}  & **GRAND+** & **81.931\({}^{*}_{*,03}\)** & **83.45\({}^{*}_{*,01}\)** & **81.259\({}^{*}_{*,02}\)** & **82.47\({}^{*}_{*,12}\)** & **83.15\({}^{*}_{*,01}\)** & **81.52\({}^{*}_{*,01}\)** & **81.09\({}^{*}_{*,01}\)** \\ \cline{2-10}  & **ResGCN** & 76.30 & 77.30 & 76.20 & 77.60 & 73.30 & 31.90 & 31.90 \\  & **ResGCN+** & **77.801\({}^{*}_{*,007}\)** & **78.701\({}^{*}_{*,007}\)** & **78.601\({}^{*}_{*,007}\)** & **76.903\({}^{*}_{*,007}\)** & **76.801\({}^{*}_{*,01}\)** & **73.60\({}^{*}_{*,01}\)** & **33.60\({}^{*}_{*,01}\)** \\ \cline{2-10}  & **DenseGCN** & 76.60 & 78.50 & 76.00 & - & - & - & - \\  & **DenseGCN+** & **78.701\({}^{*}_{*,104}\)** & **78.700\({}_{*,020}\)** & **76.901\({}^{*}_{*,104}\)** & - & - & - & - \\ \cline{2-10}  & **GCNII** & 76.40 & 81.90 & 81.50 & 84.80\({}^{*}_{*,01}\)** & 84.703\({}^{*}_{*,01}\)** & **85.20\({}^{*}_{*,01}\)** & **85.040\({}^{*}_{*,001}\)** & **86.30\({}^{*}_{*,007}\)** & **85.60\({}^{*}_{*,007}\)** \\ \cline{2-10}  & **GCNII+** & **84.706\({}^{*}_{*,007}\)** & **84.802\({}^{*}_{*,007}\)** & **84.703\({}^{*}_{*,01}\)** & **85.20\({}^{*}_{*,01}\)** & **85.20\({}^{*}_{*,01}\)** & **85.040\({}^{*}_{*,001}\)** & **85.60\({}^{*}_{*,007}\)** & **85.60\({}^{*}_{*,007}\)** \\ \hline \multirow{4}{*}{**Citeseer**} & GCN & 70.20 & 62.50 & 62.90 & 21.00 & 11.90 & 22.90 & 19.80 \\  & **GCN+** & **72.904\({}^{*}_{*,207}\)** & **67.300\({}^{*}_{*,207}\)** & **72.000\({}^{*}_{*,01}\)** & **54.703\({}^{*}_{*,31}\)** & **70.303\({}^{*}_{*,24}\)** & **48.42\({}^{*}_{*,25}\)** & **46.60\({}^{*}_{*,28}\)** \\ \cline{2-10}  & **GCN+** & **72.905\({}^{*}_{*,207}\)** & **67.300\({}^{*}_{*,207}\)** & **72.000\({}^{*}_{*,01}\)** & **54.703\({}^{*}_{*,31}\)** & **70.303\({}^{*}_{*,24}\)** & **48.42\({}^{*}_{*,25}\)** & **46.60\({}^{*}_{*,28}\)** \\ \hline \multirow{4}{*}{**Citeseer**} & GAT & 71.70 & 58.60 & 26.60 & 18.10 & - & - & - \\  & **GAT+** & **73.004\({}^{*}_{*,107}\)** & **69.501\({}^{*}_{*,017}\)** & **47.606\({}^{*}_{*,01}\)** & **31.801\({}^{*}_{*,17}\)** & **31.303\({}^{*}_{*,17}\)** & **30.605\({}^{*}_{*,01}\)** & **29.30\({}^{*}_{*,01}\)** \\ \cline{2-10}  & **GRAND+** & 71.94 & 72.58 & 73.87 & 75.00 & 75.16 & 72.90 & 69.52 \\  & **GRAND+** & **72.262\({}^{*}_{*,027}\)** & **73.550\({}^{*}_{*,017}\)** & **75.161\({}^{*}_{*,207}\)** & **75.650\({}^{*}_{*,05}\)** & **75.52\({}^{*}_{*,007}\)** & **74.521\({}^{*}_{*,027}\)** & **72.26\({}^{*}_{*,74}\)** \\ \cline{2-10}  & **ResGCN** & 67.10 & 66.00 & 63.60 & 65.00 & 62.50 & 18.80 & 80.00 & 18.10 \\ \cline{2-10}  & **ResGCN+** & **68.000\({}^{*}_{*,007}\)** & **66.001\({}^{*}_{*,001}\)** & **66.000\({}^{*}_{*,02}\)** & **66.000\({}^{*}_{*,001}\)** & **66.000\({}^{*}_{*,05}\)** & **65.80\({}^{*}_{*,01}\)** & **24.30\({}^{*}_{*,207}\)** \\ \cline{2-10}  & **DenseGCN** & 67.40 & 64.00 & 62.20 & - & - & - & - & - \\
### Experimental results on graph node classification

We postulate that by mitigating the over-smoothing issue, we can leverage the depth of GNN models to effectively capture complex feature representations in graph data. As shown in Table 1, we investigate the graph node classification accuracy as we increase the number of GNN layers by six benchmark GNN models and their corresponding plug-in models (indicated by '+' at the end of each GNN model name) with the DC-layer. The results demonstrate that: (1) the new GNN models generated from the _GNN-PDE-COV_ framework have achieved SOTA in _Cora_ (86.30% by GCNII+), _Citeseer_ (75.65% by GRAND+), and _Pubmed_ (80.10 % by GCNII+); (2) all of new GNN models outperforms their original counterparts with significant improvement in accuracy; (3) the new GNN models exhibit less sensitivity to the increase of model depth compared to current GNN models; (4) the new GNN models are also effective in resolving the gradient explosion problem [30] (e.g, the gradient explosion occurs when training GAT on all involved datasets with deeper than 16 layers, while our GAT+ can maintain reasonable learning performance even reaching 128 layers.)

It is important to note that due to the nature of the graph diffusion process, graph embeddings from all GNN models (including ours) will eventually become identical after a sufficiently large number of layers [11]. However, the selective diffusion mechanism (i.e., penalizing excessive diffusion across communities) provided by our _GNN-PDE-COV_ framework allows us to control the diffusion patterns and optimize them for specific graph learning applications.

### Application for uncovering the propagation mechanism of pathological events in AD

_First_, we evaluate the prediction accuracy between the ground truth and the estimated concentration values by our _X-FlowNet_ and six benchmark GNN methods. The statistics of MAE (mean absolute error) by _X-FlowNet_, GCN, GAT, GRAND, ResGCN, DenseGCN and GCNII, at different noise levels on the observed concentration levels, are shown in Fig. 5 (a). It is clear that our _X-FlowNet_ consistently outperforms the other GCN-based models in all imaging modalities.

_Second_, we have evaluated the potential of disease risk prediction and presented the results in Table 4 in Supplementary document, where our _GNN-PDE-COV_ model not only achieved the highest diagnostic accuracy but also demonstrated a significant improvement (paired _t_-test \(p<0.001\)) in disease risk prediction compared to other methods. These results suggest that our approach holds great clinical value for disease early diagnosis.

_Third_, we examine the spreading flows of tau aggregates in CN (cognitively normal) and AD groups. As the inward and outward flows shown in Fig. 5(b), it is evident that there are significantly larger amount of tau spreading between sub-cortical regions and entorhinal cortex in CN (early sign of AD onset) while the volume of subcortical-entorhinal tau spreading is greatly reduced in the late stage of AD. This is consistent with current clinical findings that tau pathology starts from sub-cortical regions and then switches to cortical-cortical propagation as disease progresses [28]. However, our _Tau-FlowNet_ offers a fine-granularity brain mapping of region-to-region spreading flows over time, which provides a new window to understand the tau propagation mechanism in AD etiology [14].

## 4 Conclusion

In this work, we present the _GNN-PDE-COV_ framework to re-think and re-design GNN models with great mathematical insight. On top of this, we devise the selective inductive bias to address the over-smoothing problem in GNN and develop new GNN model to predict the pathology flows _in-vivo_ via longitudinal neuroimages. Future work may involve exploring innovative graph regularization techniques and conducting further validation on a broader range of graph-based learning tasks.
## 5 Solving variational problems: From objective functional to E-L equations

### Step-by-step derivation of min-max optimization in Section 2.2.1

By substituting Eq. 2 into Eq. 1 in the main manuscript, we can obtain the objective function of subscript \(z\) (we temporarily drop \(i\) for clarity):

\[\mathcal{J}(z) =\max_{|z|\leq\mathbf{1}}\left\|\frac{\lambda}{2}z\nabla_{\mathcal{ G}}x\right\|_{2}^{2}+\lambda z\nabla_{\mathcal{G}}(x^{0}-\frac{\lambda}{2}z \nabla\mathcal{G}x) \tag{8}\] \[=\max_{|z|\leq\mathbf{1}}-\frac{\lambda^{2}}{4}z\nabla_{ \mathcal{G}}xz\nabla_{\mathcal{G}}x+\lambda z\nabla_{\mathcal{G}}x^{0} \tag{9}\]

Next, we convert Eq. 9 into a minimization problem as follows:

\[z=\begin{array}{c}\operatorname*{arg\,min}\\ |z|\leq\mathbf{1}\end{array}\quad z\nabla_{\mathcal{G}}xz\nabla_{\mathcal{G}}x -\tfrac{4}{\lambda}z\nabla_{\mathcal{G}}x^{0} \tag{10}\]

By letting the derivative with respect to \(z_{i}\) to zero, we have the following equation

\[\nabla_{\mathcal{G}}xz\nabla_{\mathcal{G}}x=\frac{4}{\lambda}\nabla_{ \mathcal{G}}x^{0} \tag{11}\]

Since \(z\) might be in high dimensional space, solving such a large system of linear equations under the constraint \(|z|\leq 1\) is oftentimes computationally challenging. In order to find a practical solution for \(z\) that satisfies the constrained minimization problem in Eq. 10, we resort to the majorization-minimization (MM) method [16]. First, we define:

\[M(z)=z\nabla_{\mathcal{G}}xz\nabla_{\mathcal{G}}x-\frac{4}{\lambda}z\nabla_{ \mathcal{G}}x^{0} \tag{12}\]

By setting \(z^{l}\) as point of coincidence, we can find a separable majorizer of \(M(z)\) by adding the non-negative function

\[(z-z^{l})^{\intercal}(\beta I-\nabla_{\mathcal{G}}x\nabla_{ \mathcal{G}}x^{\intercal})(z-z^{l}) \tag{13}\]

to \(M(z)\), where \(\beta\) is greater than or equal to the maximum eigenvalue of \(\nabla_{\mathcal{G}}x\nabla_{\mathcal{G}}x^{\intercal}\). Note, to unify the format, we use the matrix transpose property in Eq. 13. Therefore, a majorizer of \(M(z)\) is given

Figure 5: (a) Prediction accuracy by _X-FlowNet_ and six benchmark GNN models w.r.t. various noise levels. (b) The subcortical\(-\)cortical tau flows are profound in CN. But in AD, there is a diminished extent of such flows.
by:

\[M(z)+(z-z^{l})^{\intercal}(\beta I-\nabla_{\mathcal{G}}x\nabla_{\mathcal{G}}x^{ \intercal})(z-z^{l}) \tag{14}\]

And, using the MM approach, we can obtain the update equation for \(z\) as follows:

\[\begin{split} z^{l+1}&=\operatorname*{arg\,min}_{|z |\leq 1}(M(z)+(z-z^{l})^{\intercal}(\beta I-\nabla_{\mathcal{G}}x\nabla \mathcal{G}x^{\intercal})(z-z^{l}))\\ &=\operatorname*{arg\,min}_{|z|\leq 1}(\beta z^{\intercal}z-2( \nabla_{\mathcal{G}}(\frac{2}{\lambda}x^{0}-\nabla\mathcal{G}xz^{l})+\beta z ^{l})^{\intercal}z)\\ &=\operatorname*{arg\,min}_{|z|\leq 1}(z^{\intercal}z-2(\frac{1}{ \beta}\nabla_{\mathcal{G}}(\frac{2}{\lambda}x^{0}-\nabla\mathcal{G}xz^{l})+z ^{l})^{\intercal}z)\\ &=\operatorname*{arg\,min}_{|z|\leq 1}(z^{\intercal}z-2b^{\intercal}z) \end{split} \tag{15}\]

where \(b=z^{l}+\frac{1}{\beta}\nabla_{\mathcal{G}}(\frac{2}{\lambda}x^{0}-\nabla_{ \mathcal{G}}xz^{l})\).

Then, the next step is to find \(z\in\mathcal{R}^{N}\) that minimizes \(z^{\intercal}z-2bz\) subject to the constraint \(|z|\leq 1\). Let's first consider the simplest case where \(z\) is a scalar:

\[\operatorname*{arg\,min}_{|z|\leq 1}\quad z^{2}-2bz \tag{16}\]

The minimum of \(z^{2}-2bz\) is at \(z=b\). If \(b\leq 1\), then the solution is \(z=b\). If \(|b|\geq 1\), then the solution is \(z=sign(b)\). We can define the clipping function as:

\[clip(b,1):=\left\{\begin{array}{cc}&b\quad|b|\leq 1\\ &sign(b)\quad|b|\geq 1\end{array}\right. \tag{17}\]

as illustrated in the middle of Fig. 3 in the main text, then we can write the solution to Eq. 16 as \(z=clip(b,1)\).

Note that the vector case Eq. 15 is separable - the elements of \(z\) are uncoupled so the constrained minimization can be performed element-wise. Therefore, an update equation for \(z\) is given by:

\[z^{l+1}=clip(z^{l}+\frac{1}{\beta}\nabla_{\mathcal{G}}(\frac{2}{\lambda}x^{0} -\nabla_{\mathcal{G}}xz^{l}),1) \tag{18}\]

where \(l\) denotes the index of the network layer, the representation of \((l+1)^{th}\) is given by Eq. (1) in the main manuscript. Because the optimization problem is convex, the iteration will converge from any initialization. We may choose, say \(z^{0}=0\). We call this the iterative _diffusion-clip (DC) algorithm_.

This algorithm can also be written as

\[\begin{split} x^{l+1}&=x^{0}-\frac{\lambda}{2}\nabla _{\mathcal{G}}{}^{\intercal}z^{l}\\ z^{l+1}&=clip\left(z^{l}+\frac{2}{\beta\lambda} \nabla_{\mathcal{G}}x^{l+1},1\right).\end{split} \tag{19}\]

By scaling \(z\) with a factor of \(\lambda/2\), we have the following equivalent formulations:

\[\begin{split} x^{l+1}&=x^{0}-\nabla_{\mathcal{G}}{} ^{\intercal}z^{l}\\ z^{l+1}&=clip\left(z^{(i)}+\frac{1}{\beta}\nabla_{ \mathcal{G}}x^{l+1},\frac{\lambda}{2}\right)\end{split} \tag{20}\]

We summarize the process of the diffusion-clip (DC) layer in Algorithm 1 (it is similar to the iterative shrinkage threshold algorithm [12]):
The objective function:

\[\mathcal{J}(x)=\min_{x}\bigl{(}\bigl{\|}x-x^{(0)}\bigr{\|}_{2}^{2}+\mathcal{J}_{ TV}(x)\bigr{)}\]

can be minimized by alternating the following two steps:

\[x^{l} =x^{0}-\nabla_{\mathcal{G}}x^{\intercal}z^{l-1}\] \[z^{l} =clip\left(z^{l-1}+\tfrac{1}{\beta}\nabla_{\mathcal{G}}x^{l}, \tfrac{\lambda}{2}\right)=clip\left(z^{l-1}+\tfrac{2}{\beta\lambda}\nabla_{ \mathcal{G}}x^{l},1\right)\]

for \(l\geq 1\) with \(z^{0}=0\) and \(\beta\geq maxeig(\nabla_{\mathcal{G}}x^{\intercal}\nabla_{\mathcal{G}}x)\)

### The step-by-step derivation of min-max optimization schema in Section 2.2.2

According to the introduction of Section 2.2.2 (Eq. 4 and Eq. 5) in the main manuscript, we summarize the following equations,

\[\left\{\begin{array}{c}\frac{dx}{dt}+div(q)=0\\ u_{i}=\phi(x_{i})\\ q=\alpha\otimes\nabla u\\ \Delta u=div(\nabla u)\end{array}\right.\xrightarrow{\text{ derive}}\left\{\begin{array}{c}\frac{dx}{dt}=-div(q)\\ \frac{du}{dt}=-\phi^{-1}div(q)\\ \frac{du}{dt}=-\phi^{-1}div(\alpha\otimes q)\\ \frac{du}{dt}=-\phi^{-1}(\alpha\otimes\Delta u)\end{array}\right. \tag{21}\]

Since the PDE in Eq. 5 in the main manuscript is equivalent to the E-L equation of the quadratic functional \(\mathcal{J}(u)=\min\limits_{u}\int_{\mathcal{G}}\alpha\otimes|\nabla_{ \mathcal{G}}u|^{2}du\) (after taking \(\phi\) away), we propose to replace the \(\ell_{2}\)-norm integral functional \(\mathcal{J}(u)\) with TV-based counterpart

\[\mathcal{J}_{TV}(u)=\min\limits_{u}\int_{\mathcal{G}}\alpha\otimes|\nabla_{ \mathcal{G}}u|du \tag{22}\]

We then introduce an auxiliary matrix \(f\) to lift the undifferentiable barrier, and reformulate the TV-based functional as a dual min-max functional

\[\mathcal{J}_{TV}(u,f)=\min\limits_{u}\max\limits_{f}\int_{\mathcal{G}}\alpha \otimes f(\nabla_{\mathcal{G}}u)du \tag{23}\]

where we maximize \(f\) such that \(\mathcal{J}_{TV}(u,f)\) is close enough to \(\mathcal{J}_{TV}(u)\). Using Gateaux variations, we assume \(u\to u+\varepsilon a\), \(f\to f+\varepsilon b\), and the directional derivatives in the directions \(a\) and \(b\) defined as \(\frac{d\mathcal{J}}{d\varepsilon}(u+\varepsilon a)\big{|}_{\varepsilon\to 0}\) and \(\frac{d\mathcal{J}}{d\varepsilon}(f+\varepsilon b)\big{|}_{\varepsilon\to 0}\). Given a functional \(\mathcal{J}_{TV}(u,f)\), its Gateaux variations is formulated by:

\[\mathcal{J}_{TV}(u+\varepsilon a,f+\varepsilon b)=\int\alpha \otimes[(f+\varepsilon b)\cdot(\nabla u+\varepsilon\nabla a)]du \tag{24}\] \[\Rightarrow\frac{\partial\mathcal{J}}{\partial\varepsilon}\bigg{|} _{\varepsilon\to 0}=\int\alpha\otimes[(f\cdot\nabla a)+(\nabla ub)]\ du\] \[\Rightarrow\frac{\partial\mathcal{J}}{\partial\varepsilon}\bigg{|} _{\varepsilon\to 0}=\alpha\otimes f\cdot a-\int\alpha\otimes(a\cdot\nabla f)du+ \int\alpha\otimes(b\nabla u)du\]

Since we assume either \(u\) is given at the boundary (Dirichlet boundary condition), the boundary term \(\alpha\otimes f\cdot a\) can be dropped. After that, the derivative of \(\mathcal{J}_{TV}(u,f)\) becomes:

\[\frac{\partial\mathcal{J}}{\partial\varepsilon}\bigg{|}_{\varepsilon\to 0}=-\int \alpha\otimes(\nabla f\cdot a+\nabla u\cdot b) \tag{25}\]
Since the dummy functional \(a\) and \(b\) are related to \(u\) and \(f\) respectively, the E-L equation from the Gateaux variations in Eq. 25 leads to two coupled PDEs:

\[\left\{\begin{array}{l}\max\limits_{f}\frac{df}{dt}=\alpha\otimes\nabla_{\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
* Following the region parcellations, we calculate the regional concentration level (the Cerebellum as the reference) of the amyloid, FDG, CoTh and tau pathologies for each brain region (red arrows in Fig. 6), yielding the input \(x\in\mathcal{R}^{160}\) for training _X-FlowNet_, respectively.

Following the clinical outcomes, we partition the subjects into the cognitive normal (CN), early-stage mild cognitive impairment (EMCI), late-stage mild cognitive impairment (LMCI), and AD groups. To facilitate population counts, we regard CN and EMCI as "CN-like" group, while LMCI and AD as "AD-like" groups. Table 3 summarizes the statistics of the two datasets.

### Experiments on node classification

Fig 7 presents the performance of different evaluation criteria (accuracy, precision, and F1-score) across different network layers for node classification by benchmark GNN model (patterned in dash lines) and the counterpart novel GNN model from our _GNN-PDE-COV_ framework (patterned by solid lines), where each row is associated with a specific instance of GNN model. It is evident that our proposed _GNN-PDE-COV_ consistently outperforms other methods across different layers, with significantly enhanced degrees in accuracy, precision, and F1-score. Moreover, the GNN model yielded from our _GNN-PDE-COV_ framework consistently achieves the highest accuracy on all three datasets. Overall, these results demonstrate the state-of-the-art performance by our _GNN-PDE-COV_ framework in graph node classification.

The effect of anti-smoothing by clip operation is shown in Fig. 8. To set up the stage, we put the spotlight on the links that connect two nodes with different categorical labels. In this context, 2,006

\begin{table}
\begin{tabular}{l c c c c c c} \hline \hline Algorithm & Optimizer & Learning rate & Weight decay & Hidden layer & Dropout & Epoch \\ \hline _GCN_ & Adam & 0.01 & \(5\times 10^{-4}\) & \(M\to 2^{i}\rightarrow...\to 2^{4}\to C\) & 0.5 & 1500 \\ \hline _GAT_ & Adam & 0.001 & \(5\times 10^{-4}\) & head=8, \(M\to 2^{i}...\to C\) & 0.6 & 2000 \\ \hline _RGCN_ & Adam & 0.005 & \(5\times 10^{-4}\) & hidden dimension=64 & 0.1 & 2500 \\ \hline _DGCN_ & Adam & 0.001 & \(5\times 10^{-4}\) & hidden dimension=64 & 0.1 & 2500 \\ \hline _GRAND_ & Adam & 0.01 & \(5\times 10^{-4}\) & hidden dimension=16 & 0.5 & 200 \\ \hline _GCNII_ & Adam & 0.005 & \(5\times 10^{-4}\) & hidden dimension=128 & 0.6 & 2000 \\ \hline \hline _GCN_ & Adam & 0.001 & \(5\times 10^{-4}\) & hidden dimension=16 & 0.2 & 500 \\ \hline _GAT_ & Adam & 0.001 & \(5\times 10^{-4}\) & head=8, hidden dimension=4 & 0.5 & 800 \\ \hline _RGCN_ & Adam & 0.001 & \(5\times 10^{-4}\) & hidden dimension=16 & 0.1 & 500 \\ \hline _DGCN_ & Adam & 0.01 & \(5\times 10^{-4}\) & hidden dimension=8 & 0.1 & 500 \\ \hline _GCNII_ & Adam & 0.005 & \(5\times 10^{-4}\) & hidden dimension=16 & 0.6 & 1500 \\ \hline _GRAND_ & Adam & 0.01 & \(5\times 10^{-4}\) & hidden dimension=16 & 0.5 & 500 \\ \hline _X-FlowNet_ & Adam & 1e-4/3e-3 & \(1\times 10^{-5}\) & hidden dimension=16 & 0.5 & 500 \\ \hline _PDENet_ & Adam & 0.01 & \(1\times 10^{-5}\) & hidden dimension=16 & 0.5 & 500 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Parameters setting on Citation network (top) and ADNI data (bottom). \(M\) denotes the feature dimension and \(C\) denotes the number of classes. For _Cora_ dataset, we set \(i=4\) when network layer \(L=2\), \(i=8\) if \(L=4\), \(i=10\) if \(L=8,16,32,64,128\). For _Citeseer_ dataset, we set \(i=4\) when network layer \(L=2\), \(i=8\) if \(L=4\), \(i=11\) if \(L=8,16,32,64,128\). For Pubmed dataset, we set \(i=4\) when network layer \(L=2\), \(i=8\) if \(L=4,8,16,32,64,128\). The hidden dimension of \(l^{th}\) is twice that of layer \((l-1)^{th}\). Take Cora as an example (8 layers), the dimension of the hidden layer is: 1433 \(\rightarrow\) 1024 \(\rightarrow\)512 \(\rightarrow\) 256\(\rightarrow\)128\(\rightarrow\)64\(\rightarrow\)32\(\rightarrow\)16\(\rightarrow\)7. After exceeding 8 layers, the number of hidden layers is doubled according to the total network layer.

\begin{table}
\begin{tabular}{l c c c|l c c} \hline \hline \multicolumn{4}{c|}{**Node classification (Citation)**} & \multicolumn{3}{c}{**Application on flow prediction (ADNI)**} \\ \hline \multirow{2}{*}{Dataset} & \multicolumn{3}{c}{Description} & \multicolumn{1}{c}{Features} & \multicolumn{1}{c}{\# of subjects (CN/AD)} \\ \cline{2-7}  & Classes & Nodes & Edges & Features & Amyloid (160) & 304/83 \\ \hline _Cora_ & 7 & 2708 & 5429 & 1433 & Tau (160) & 124/37 \\ _Citeseer_ & 6 & 3327 & 4732 & 3703 & FDG (160) & 211/63 \\ _Pubmed_ & 3 & 19717 & 44338 & 500 & Cortical thickness (160) & 359/110 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Dataset statistics.
links from _Cora_, 2,408 links from _Citeseer_, and 17,518 links from _Pubmed_ datasets are selected, called inter-class links. For each inter-class link, we calculate node-to-node similarity in terms of Pearson's correlation between two associated graph embedding vectors 4 by benchmark methods (in red) and the counterpart GNN models derived from _GNN-PDE-COV_ framework (in green). We find that (1) more than 70% nodes are actually associated with inter-class links which confirms the hypothesis of over-smoothing in Fig. 1 of our manuscript; (2) Our novel GNN models have the ability to learn feature representations that better preserve the discriminative power for node classification (as indicated by the distribution of node-to-node similarity shifting towards the sign of anti-correlation).

Footnote 4: the learned feature representations for node classification

### Application on uncovering the propagation mechanism of pathological events in AD

_Firstly_, we examine the prediction accuracy for each modality of concentration (tau, amyloid, FDG, CoTh) level at different noise levels. Specifically, to evaluate the robustness of our _X-FlowNet_ model to noise, we conducted an experiment by adding uncorrelated additive Gaussian noise levels with standard deviation ranging from 0.02 to 1 to the observed concentration levels of tau, amyloid, FDG, and CoTh. We then evaluated the prediction accuracy (MAE) using 5-fold cross-validation. The prediction results, as shown in Fig. 9, indicate that our _X-FlowNet_ model is less sensitive to noise added to the imaging features than all other counterpart GNN methods.

_Secondly_, we conduct an ablation study to compare our _X-FlowNet_ model with _PDENet_ (marked as #7 in Fig. 9). Our model, which is in a GAN architecture and incorporates a TV constraint to avoid over-smoothing, integrates the two steps of estimating the PEF and uncovering the spreading flows into a unified neural network, resulting in significantly improved prediction accuracy compared to _PDENet_.

_Thirdly_, we perform a disease risk prediction experiment, which can be regarded as a graph classification problem. We assume that we have baseline amyloid, tau, FDG, and CoTh scans, and evaluate the prediction accuracy, precision and F1-score of various models in forecasting the risk of

Figure 6: General workflows for processing T1-weighted image (yellow arrows), diffusion-weighted image (white arrows), and PET images (red arrows). The output is shown at the bottom right, including the brain network, and regional concentration level of amyloid, FDG, CoTh and tau aggregates.
developing AD. We consider two dichotomous cases: one included only AD vs. CN groups and the other involved AD/LMCI vs. CN/EMCI. The results of the mean of 5-fold cross-validation are shown in Table 4. Our _GNN-PDE-COV_ outperforms all other methods in terms of accuracy, precision and F1-score indicated by an asterisk (**'*') at the significance level of 0.001.

_Fourthly_, we examine the propagation pattern of tau spreading flows on an individual basis (Fig. 10). _First_, we visualize the top flows (ranked in terms of flow volume) uncovered in a CN subject (Fig. 10(a)). It is apparent that subcortex-cortex flows are the predominant patterns, where most of the tau aggregates spread from subcortical regions (globus pallidus, hippocampus, and putamen) to the temporal lobe, limbic lobe, parietal lobe, and insula lobe. Note, we find inferior temporal gyrus (\(t_{6}\)) and entorhinal cortex (\(t_{8}\)) are actively involved in the subcortex-cortex flows, which are the footprints of early stage tau propagation frequently reported in many pathology studies [28; 40]. _Second_, we visualize the top flows uncovered in an AD subject (Fig. 10(b)). It is apparent that the propagation of tau is restricted on the brain cortex, mainly spreading from temporal lobe regions to other regions (such as frontal lobe, limbic lobe and occipital lobe), which is aligned with current clinical and pathology findings that predominant amount of tau aggregates propagate throughout brain cortex in the late stage of AD.

### Discussion and limitations

_Discussion._ In our experiments, we found adding DC layer right after every FC layer usually does not yield best performance. Instead, we empirically set to add DC layer from the first several FC layers. For example, we add DC layer after the \(3^{rd}\) FC layer in an 8-layer GNN model, after the \(5^{th}\) FC layer in a 16-layer GNN model, and after \(8^{th}\) FC layer in a GNN model with more than 16 layers. One possible explanation is that the clip operation in DC layer depends on a good estimation of cap \(b\) in Eq. 3 (in the main manuscript). Given that the estimation of \(b\) may lack stability during the initial stages of graph learning, it can be advantageous to postpone the clip operation from an engineering perspective. However, delaying the addition of the DC layer too much can result in missed opportunities to address the problem of over-smoothing.

Regarding the computational time, we record the additional computational time of training our DC layer on different datasets. Specifically, the extra training time is 2.2 ms/epoch in _Cora_, 9.8

\begin{table}
\begin{tabular}{c|c|c c c c|c c c c|c c} \hline \hline
**Tau** & Unit (\%) & GCN & **GCN+** & GAT & **GAT+** & GCNII & **GCNI+** & RGCN & **RGCN+** & DGGCN+** & GRAND & **GRAND+** \\ \hline AD/LMCI & Precision & 80.15 & **90.03**(*) & 69.91 & **86.18*** & 83.93 & **90.03**(*) & 84.64 & **89.46*** & 84.03 & **91.58*** & 87.95 & **88.22*** \\ vs. & Accuracy & 82.30 & **88.74*** & 81.05 & **87.50*** & 83.79 & **88.75*** & 86.03 & **90.00*** & 85.54 & **91.25*** & 88.75 & **86.12*** \\ CN/EMCI & F1-score & 75.55 & **84.49*** & 72.87 & **84.72*** & 88.44*** & 83.18 & **85.54*** & 82.45 & **91.39*** & 88.41 & **89.44*** \\ \hline AD & Precision & 89.29 & **91.92*** & 97.82 & 96.90*** & 83.65 & **88.52*** & 92.61 & **95.72*** & 92.61 & **95.91*** & 91.77 & **95.76*** \\ vs. & Accuracy & 86.64 & **90.91*** & 84.86 & **88.14*** & 76.84 & 86.36 & 91.01 & 97.45*** & 91.97 & **95.65*** & 90.91 & **95.45*** \\ CN & F1-score & 86.54 & **90.26*** & 83.39 & **89.71*** & **76.1*** & 75.11 & **84.68*** & 90.45 & 95.32*** & 90.54 & 95.55*** & 88.86 & **95.38*** \\ \hline
**Amyloid** & Unit (\%) & GCN & **GCN+** & GAT & **GAT+** & GCNII & **GCNI+** & RGCN & **RGCN & **DGCN+** & GRAND & **GRAND+** \\ \hline AD/LMCI & Precision & 76.36 & **83.78*** & 67.73 & **71.79*** & 60.87 & **0.01**(*) & 72.53 & **83.21*** & 74.92 & **60.17*** & 79.00 & **79.93*** \\ vs. & Accuracy & 76.40 & **79.44*** & 75.43 & **75.37*** & **75.43*** & 74.61*** & 75.90 & **78.50*** & 76.92 & 75.77 & **80.37** & **81.31*** \\ CN/EMCI & F1-score & 70.33 & **72.58*** & 67.66 & **66.93*** & 63.57 & **67.31*** & 70.66 & **70.67*** & **70.62*** & **76.62*** & **77.79*** & **79.63*** \\ \hline AD & Precision & 81.58 & **88.37*** & 81.54 & **87.98*** & 70.59 & **79.98*** & 85.75 & **93.09*** & 83.87 & **90.56*** & 65.53 & **89.62*** \\ vs. & Accuracy & 80.77 & **88.10*** & 80.78 & **85.10*** & 75.02 & **81.24*** & 85.56 & **92.86*** & 85.29 & **90.48*** & 80.95 & **87.80*** \\ CN & F1-score & 78.14 & **87.68*** & 78.07 & **87.98*** & 65.87 & **77.42*** & 85.34 & **92.92*** & **82.30** & **92.07*** & 72.43 & **88.22*** \\ \hline
**FDG** & Unit (\%) & GCN & GCN & **GAT+** & GCNII & **GCNI+** & RGCN & **RGCN & DGCN+** & GRAND & **GRAND+** \\ \hline AD/LMCI & Precision & 84.83 & **67.20*** & 55.86 & **62.90*** & 60.88 & **93.04*** & 90.45 & **85.14*** & 90.45 & **85.14*** & 91.38 & **86.25*** \\ vs. & Accuracy & 73.17 & **76.00*** & 72.17 & **77.00*** & 71.78 & **74.54*** & 70.98 & **74.26*** & 70.98 & **74.26*** & 71.10 & **75.00*** \\ CN/EMCI & F1-score & 89.46 & **81.56*** & 69.12 & **65.90*** & 61.02 & **69.07*** & 85.89 & **63.29*** & 85.86 & **63.29*** & 59.42 & **64.29*** \\ \hline AD & Precision & 81.31 & **87.25*** & 61.90 & **82.33*** & 74.11 & **81.06*** & 90.57 & **89.77** & 89.50*** & 63.84 & **81.77*** & 70.91 & **72.24****** \\ vs. & Accuracy & 82.17 & **84.62*** & 72.82 & **78.95*** & 79.55 & **82.05*** & 73.35 & **79.58*** & 73.80*** & **80.11*** & 84.21 & **86.32*** \\ CN & F1-score & 79.40 & **82.04*** & 64.23 & **69.66*** & 73.38 & **88.09*** & 62.77 & **75.98*** & **63.92*** & **76.58*** & **76.59*** & **78.06*** \\ \hline \hline CoIn & Unit (\%) & GCN & **GCN-** & GAT & **GAT+** & GCNII & **GCNI+** & RGCN & **RGCN+** & GRAND & **GRAND+** \\ AD/LMCI & Precision & 74.85 & **74.74*** & 62.63 & **67.15*** & 62.63 & **74.71*** & 62.63 & **68.77*** & 62.63 & **64.59****** & 63.81 & **68.77*** \\ vs. & Accuracy & 80.68 & **82.32*** & 9.10 & **79.19*** & **79.34*** & 9.10 & **80.37*** & 99.10 & **80.37*** & **97.88** & **82.93*** \\ CN/EMCI & F1-score & 73.55 & **75.93*** & 69.89 & **70.44*** & 69.89 & **72.72*** & 69.89 & **75.19*** & 69.89 & **71.62*** & 70.94 & **75.19*** \\ \hline AD & Precision & 83.45 & **85.77*** & 71.24 & **72.80*** & **78.0*** & 76.14 & **79.84*** & 65.04 & **
ms/epoch in _Citeseer_, 7.8 ms/epoch in _Pubmed_, and 0.3 ms/epoch in _ADNI_, respectively, where the data descriptions are listed in Table 3. It is apparent that the TV-based constraint effectively addresses the over-smoothing issue in GNN without imposing a significant computational burden.

_Limitations._ Our current graph learning experiments are limited to citation networks. In the future, we will evaluate our _GNN-PDE-COV_ framework on other graph datasets such as drug medicine and protein networks.

_Societal impact._ Our major contribution to the machine learning field is a novel research framework which allows us to develop new GNN models with a system-level understanding. We have provided a new approach to address the common issue of over-smoothing in GNN with a mathematical guarantee. From the application perspective, the new deep model for uncovering the _in-vivo_ propagation flows has great potential to establish new underpinning of disease progression and disentangle the heterogeneity of diverse neurodegeneration trajectories.
