Feature Transportation Improves Graph Neural Networks
Graph neural networks (GNNs) have shown remarkable success in learning
representations for graph-structured data. However, GNNs still face challenges
in modeling complex phenomena that involve feature transportation. In this
paper, we propose a novel GNN architecture inspired by
Advection-Diffusion-Reaction systems, called ADR-GNN. Advection models feature
transportation, while diffusion captures the local smoothing of features, and
reaction represents the non-linear transformation between feature channels. We
provide an analysis of the qualitative behavior of ADR-GNN, that shows the
benefit of combining advection, diffusion, and reaction. To demonstrate its
efficacy, we evaluate ADR-GNN on real-world node classification and
spatio-temporal datasets, and show that it improves or offers competitive
performance compared to state-of-the-art networks.
Moshe Eliasof, Eldad Haber, Eran Treister
# ADR-GNN: Advection-Diffusion-Reaction Graph Neural Networks

###### Abstract

Graph neural networks (GNNs) have shown remarkable success in learning representations for graph-structured data. However, GNNs still face challenges in modeling complex phenomena that involve advection. In this paper, we propose a novel GNN architecture based on Advection-Diffusion-Reaction systems, called ADR-GNN. Advection models the directed transportation of information, diffusion captures the local smoothing of information, and reaction represents the non-linear transformation of information in channels. We provide an analysis of the qualitative behavior of ADR-GNN, that shows the benefit of combining advection, diffusion, and reaction. To demonstrate its efficacy, we evaluate ADR-GNN on real-world node classification and spatio-temporal datasets, and show that it improves or offers competitive performance compared to state-of-the-art networks.

## 1 Introduction

Recently, GNNs have been linked to ordinary and partial differential equations (ODEs and PDEs) in a series of works [103; 14; 27; 70; 90; 26; 32]. These works propose to view GNN layers as the time discretization of ODEs and PDEs, and as such they offer both theoretical and practical advantages. First, ODE and PDE based models allow to reason about the behavior of existing GNNs. For instance, as suggested in [14], it is possible to view GCN [45] and GAT [88] as discretizations of the non-linear heat equation. This observation helps to analyze and understand the oversmoothing phenomenon in GNNs [61; 62; 13]. Second, ODE and PDE based GNNs pave the path to the construction and design of GNNs that satisfy desired properties, such as energy-preservation [27; 70], attraction and repulsion forces modeling [90; 26], anti-symmetry [32], as well as reaction-diffusion systems [22]. Nonetheless, the aforementioned architectures still rely on controlled diffusion or wave propagation, as well as non-linear pointwise convolutions. Therefore, as discussed in [69], while there are methods that can alleviate oversmoothing, they may lack expressiveness. We now provide a simple example, known as the graph node feature transportation task, where diffusion, wave propagation, and reaction (pointwise) networks may fail. In this task, the goal is to gather the node information (i.e., features) from several nodes to a single node. Clearly, no diffusion process can express or model such a phenomenon, because diffusion spreads and smooths, rather than transports information [28]. Likewise, a wave-propagation approach cannot express such a phenomenon, because it lacks directionality, which is required for this task. An instance of this problem is illustrated in Figure 1, where we show the source and target node features, and possible advection weights that can achieve the desired target. Later, in Figure 2, we show that popular operators such
as diffusion or reaction cannot model the transition from the source to the target node features, while advection can. Furthermore, the concept of advection appears in many real-world problems and data, such as traffic-flow and-control [8], quantity transportation in computational biology [87], and rainfall forecasting [75]. Motivated by the previously discussed observations and examples, we propose, in addition to learning and combining diffusion and reaction terms, to develop a learnable _advection_ term that is suited to model directed information transportation from the data. The resulting architecture, called _ADR-GNN_, can therefore express various phenomena, from advection, diffusion, to pointwise reactions, as well as their compositions.

**Contributions.** The contributions of this paper are three-fold. Firstly, we develop a novel graph neural advection operator that is mass preserving, stable, and consistent with continuous advection PDEs. This operator enables the modeling of phenomena that involve the transportation of information on graphs by learning the direction of the transportation. Secondly, we propose ADR-GNN, a GNN based on learnable advection-diffusion-reaction (ADR) systems, that can express a wide range of phenomena, including learned directional information flow, diffusion, and pointwise reactions. Thirdly, we demonstrate the efficacy of ADR-GNN on node classification, and spatio-temporal forecasting datasets, achieving improved or competitive results compared to state-of-the-art models.

## 2 Related Work

**Advection-Diffusion-Reaction.** An Advection-Diffusion-Reaction system is a mathematical model that describes the simultaneous existence of three processes: (i) the advection (transport) of information in a medium, (ii) the diffusion (smoothing) of information within that medium, and (iii) pointwise (self) reactions. These systems are used to study and model a wide range of physical, chemical, and biological phenomena. For example, ADR systems can be utilized to track and estimate the location of fish swarms [2], modeling ecological trends [24], and the modeling of turbulent flames in supernovae [44]. However, the aforementioned works rely on a low-dimensional, hand-crafted, non-neural ADR system to be determined, typically by trial and error, often requiring a domain expert. In contrast, in this paper we propose to learn the ADR system for various graph types and tasks.

**Graph Neural Networks as Dynamical Systems.** Adopting the interpretation of convolutional neural networks (CNNs) as discretizations of ODEs and PDEs [72; 20; 97] to GNNs, works like GODE [103], GRAND [14], PDE-GCN\({}_{\text{D}}\)[27], GRAND++ [86] and others, propose to view GNN layers as time steps in the integration of the non-linear heat equation, allowing to control the diffusion (smoothing) in the network, to understand oversmoothing [61; 62; 13] in GNNs. Thus, works like [21; 56; 55; 26] propose to utilize a _learnable_ diffusion term, thereby alleviating oversmoothing. Other architectures like PDE-GCN\({}_{\text{M}}\)[27] and GraphCON [70] propose to mix diffusion and oscillatory processes (e.g., based on the wave equation) to avoid oversmoothing by introducing a feature energy preservation mechanism. Nonetheless, as noted in [69], besides alleviating oversmoothing, it is also important to design GNN architectures with improved expressiveness. Recent examples of such networks are [32] that propose an anti-symmetric GNN to alleviate over-squashing [3], and [90; 22] that formulate a reaction-diffusion GNN to enable non-trivial pattern growth. In this pa

Figure 1: An example of node feature transportation on a graph. Applying the advection weights in (a) to the source (b), yields the target (c). Darker edge colors in (a) indicate greater advection weights.
per, we build on the properties of ADR PDEs, that in addition to modeling diffusive and reactive processes, also allow to capture advective processes such as the transportation of node features.

On another note, CNNs and GNNs are also used to accelerate PDE solvers [66; 53; 50; 12; 73], as well as to generate [74] and compute [6] physical simulations. In this paper, we focus on the view of GNNs as the discretization of ADR PDEs, rather than using GNNs to solve PDEs.

**Advection on Graphs.** Advection is a term used in Physics to describe the transport of a substance in a medium. In the context of graphs, advection is used to express the transport of information (features) on the graph nodes. The underlying process of advection is described by a continuous PDE, and several graph discretization techniques [15; 39] are available. The advection operator has shown its effectiveness in classical (i.e., non-neural) graph methods, from blood-vessel simulations [25], to traffic flow prediction [11]. In this paper, we develop a neural advection operator that is combined with neural diffusion and reaction operators, called ADR-GNN.

## 3 Advection-Diffusion-Reaction Graph Neural Networks

In this section, we first describe the general outline of a continuous ADR system in Section 3.1, and present its graph discrete analog, named _ADR-GNN_ in Section 3.2. We then discuss each component in ADR-GNN in detail in Sections 3.3-3.4.

**Notations.** We define a graph by \(\mathcal{G}=(\mathcal{V},\mathcal{E})\), where \(\mathcal{V}\) is a set of \(n\) nodes and \(\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}\) is a set of \(m\) edges. We denote the 1-hop neighborhood of the \(i\)-th node by \(\mathcal{N}_{i}\), and the node features by \(\mathbf{U}\in\mathbb{R}^{n\times c}\), where \(c\) is the number of features. The symmetric graph Laplacian reads \(\mathbf{L}=\mathbf{D}-\mathbf{A}\), and the symmetric normalized Laplacian is given by \(\hat{\mathbf{L}}=\mathbf{D}^{-\frac{1}{2}}\mathbf{L}\mathbf{D}^{-\frac{1}{2}}\), where \(\mathbf{D}\) is the degree matrix.

### Continuous Advection-Diffusion-Reaction Systems

The continuous PDE that describes an ADR system is given by:

\[\frac{\partial U}{\partial t}=\underbrace{\nabla\cdot(VU)}_{\text{Advection}}+ \underbrace{K\Delta U}_{\text{Diffusion}}+\underbrace{f(U,X,\theta_{r})}_{ \text{Reaction}},\ X\in\Omega,\ t\in[0,T], \tag{1}\]

accompanied by initial conditions \(U(X,t=0)\) and boundary conditions. Here, \(U(X,t)=[u_{1}(X,t),\ldots,u_{c}(X,t)]:\mathbb{R}^{\Omega\times[0,T]}\to \mathbb{R}^{c}\) is a density function, written as a vector of scalar functions \(u_{s}(X,t),\ s=1,\ldots,c\), that depend on the initial location \(X\) and time \(t\). The spatial domain \(\Omega\) can be \(\mathbb{R}^{d}\) or a manifold \(\mathcal{M}\subseteq\mathbb{R}^{d}\). From a neural network perspective, \(u_{s}\) is referred to as a channel, interacting with other channels. The left-hand side of Equation (1) is a time derivative that represents the change in features in time, as discussed in Section 2. The right-hand side includes three terms:

* **Advection.** Here, \(V\) denotes a velocity function that transports the density \(U\) in space, and \(\nabla\cdot\) is the divergence operator.
* **Diffusion.** We denote the continuous Laplacian operator by \(\Delta\). The Laplacian is scaled with a diagonal matrix \(K=\mathrm{diag}(\kappa_{1},\ldots,\kappa_{c})\in\mathbb{R}^{c\times c},\ \ \kappa_{i}\geq 0\) of non-negative diffusion coefficients, each independently applied to its corresponding channel in \(U\).
* **Reaction.** Here, \(f(U,X,\theta_{r})\) is a non-linear pointwise function parameterized by \(\theta_{r}\).

### Advection-Diffusion-Reaction on Graphs

Equation (1) is defined in the continuum. We now use a graph \(\mathcal{G}=(\mathcal{V},\mathcal{E})\) to discretize \(\Omega\). The nodes \(\mathcal{V}\) can be regarded as a discretization of \(X\), that is, the \(i\)-th node is located in \(\mathbf{X}_{i}\), and the edges \(\mathcal{E}\) represent the topology of \(\Omega\). Then, the _spatial_, graph discretization of Equation (1) is:

\[\frac{d\mathbf{U}(t)}{dt} = \mathbf{DIV}\left(\mathbf{V}\left(\mathbf{U}(t),t;\boldsymbol{ \theta}_{a}(t)\right)\mathbf{U}(t)\right)-\hat{\mathbf{L}}\mathbf{U}(t) \mathbf{K}(t;\boldsymbol{\theta}_{d}(t))+f(\mathbf{U}(t),\mathbf{X},t; \boldsymbol{\theta}_{r}(t)) \tag{2a}\] \[\mathbf{U}(0) = \mathbf{U}^{(0)}=g_{\text{in}}(\mathbf{X},\boldsymbol{\theta}_{0}) \tag{2b}\]
Here, \(\mathbf{U}(t)\in\mathbb{R}^{n\times c}\) is a matrix that describes the node features at time \(t\). The advection term depends on the velocity function \(\mathbf{V}\) parameterized by learnable weights \(\mathbf{\theta}_{a}(t)\). The precise discretization of the advection operator \(\mathbf{DIV}(\mathbf{V}\cdot)\) is discussed in Section 3.4. The diffusion is discretized using the symmetric normalized Laplacian1\(\hat{\mathbf{L}}\) that is scaled with a diagonal matrix with non-negative learnable diffusion coefficients on its diagonal \(\mathbf{K}(t;\mathbf{\theta}_{d}(t))=\mathrm{diag}(\mathrm{hardtanh}(\mathbf{\theta}_ {d}(t),0,1))\geq 0\), where \(\mathrm{hardtanh}(\mathbf{\theta}_{d}(t),0,1)\) clamps each element in \(\mathbf{\theta}_{d}\in\mathbb{R}^{c}\) to be between 0 and 1. The reaction term \(f\) from Equation (2a) is a pointwise non-linear function realized by a multilayer-perceptron (MLP) parameterized by learnable weights \(\mathbf{\theta}_{r}(t)\). To obtain initial node embedding \(\mathbf{U}^{(0)}\in\mathbb{R}^{n\times c}\) from the input features \(\mathbf{X}\in\mathbb{R}^{n\times c_{in}}\), we use a fully-connected layer \(g_{in}\) in Equation (2b).

Footnote 1: In PDE theory, the Laplacian is a negative operator, while in graph theory it is positive. Therefore it is required to multiply \(\hat{\mathbf{L}}\) by a negative sign in Equation (2a) compared to Equation (1).

In this work we focus on static and temporal node-level tasks, and we note that typically, \(c\), the number of hidden channels of \(\mathbf{U}(T)\in\mathbb{R}^{n\times c}\), is different than \(c_{out}\), the number of channels of the target output \(\mathbf{Y}\in\mathbb{R}^{n\times c_{out}}\). Therefore the output of neural network, \(\tilde{\mathbf{Y}}\), is given by

\[\tilde{\mathbf{Y}}=g_{\mathrm{out}}(\mathbf{U}(T),\mathbf{\theta}_{ \mathrm{out}})\in\mathbb{R}^{n\times c_{out}}, \tag{3}\]

where \(g_{\mathrm{out}}\) is a fully-connected layer, parameterized by learnable weights \(\mathbf{\theta}_{\mathrm{out}}\).

**The qualitative behavior of ADR-GNN.** The ADR-GNN model combines the learning of three powerful terms. Namely, the learned parameters are \(\mathbf{\theta}_{a}\) the advection parameters, \(\mathbf{\theta}_{d}\) the diffusion parameters, and \(\mathbf{\theta}_{r}\) the reaction parameters. Therefore, an ADR-GNN layer can express and model various phenomena. For example, if we set \(\mathbf{\theta}_{d}(t)=0\), then there is no diffusion in the system, and the method is dominated by advection and reaction. If on the other hand, one learns a very small advection (i.e., the learned \(\mathbf{V}\), to be discussed later, tends to retain all features in place), then a reaction-diffusion oriented system is obtained. Similarly, other combinations of advection, diffusion, and reaction can be achieved, because of the learning of the parameters of the system. Thus, ADR-GNN can be adopted to solve a host of problems, depending on dynamics and patterns mandated by the data, as we show later in our experiments in Section 4.

### From an ODE to a Graph Neural Network - Time Discretization of ADR-GNN

Equation (2) _spatially_ discretizes the PDE in Equation (1), yielding an ODE defined on the graph. The _time_ discretization of the ODE yields a sequential process that can be thought of as layers of neural networks [36; 19; 91]. That is, upon discrete time integration of Equation (2a), we replace the notion of time \(t\) with \(l\) layers, and a step size \(h\), that is a positive scalar hyperparameter.

While it is possible to use many ODE discretization methods (see, e.g., [36; 98; 19; 14] and references within), in various applications where an ADR system arises, from flow in porous media [23], to PDE-based image segmentation [89], and multiphase flow [43], an operator-splitting (OS) [4] is utilized. We therefore also use an OS time discretization for Equation (2a), that yields a graph neural ADR layer, summarized in Algorithm 1. Composing several neural ADR layers leads to ADR-GNN. We further discuss the properties of the OS approach in Appendix A. The exact discretizations of the ADR terms are derived in Section 3.4.

```
Input: Node features \(\mathbf{U}^{(l)}\in\mathbb{R}^{n\times c}\) Output: Updated node features \(\mathbf{U}^{(l+1)}\in\mathbb{R}^{n\times c}\)  Advection: \(\mathbf{U}^{(l+1/3)}=\mathbf{U}^{(l)}+h\mathbf{DIV}(\mathbf{V}(\mathbf{U}^{(l) },t;\mathbf{\theta}_{a}^{(l)})\mathbf{U}^{(l)})\)  Diffusion: \(\mathbf{U}^{(l+2/3)}=\mathrm{mat}\left(\left(\mathbf{I}+h\mathbf{K}(t;\mathbf{ \theta}_{d}^{(l)})\otimes\tilde{\mathbf{L}}\right)^{-1}\mathrm{vec}(\mathbf{U }^{(l+1/3)})\right)\)  Reaction: \(\mathbf{U}^{(l+1)}=\mathbf{U}^{(l+2/3)}+hf(\mathbf{U}^{(l+2/3)},\mathbf{U}^{( 0)},t;\mathbf{\theta}_{r}^{(l)})\)
```

**Algorithm 1** Graph Neural Advection-Diffusion-Reaction Layer

### Discretized Graph Operators

We now elaborate on the discretized graph operators utilized in our ADR-GNN, summarized in Algorithm 1. Besides the combination of the learnable advection, diffusion, and reaction terms,
which, to the best of our knowledge, was not studied in the context of GNNs, the main innovation here is the _consistent, mass preserving, and stable_ discretization of the advection operator. We further utilize recent time integration techniques of the diffusion and reaction terms, discussed below.

**Advection.** To define the graph discretized advection operator, we extend the non-learnable advection operator from [15], into a learnable, neural advection operator. Our advection operator transports node features based on learned directed edge weights (velocities) \(\{(\mathbf{V}_{i\to j},\mathbf{V}_{j\to i})\}_{(i,j)\in\mathcal{E}}\), where each \(\mathbf{V}_{i\to j},\mathbf{V}_{j\to i}\in\mathbb{R}^{c}\), such that \(0\leq\mathbf{V}_{i\to j}\leq 1\). The notation \(i\to j\) implies that the weight transfers features from the \(i\)-th to \(j\)-th node. We further demand that the outbound edge weights associated with every node, per channel, sum to 1, i.e., \(\sum_{j\in\mathcal{N}_{i}}\mathbf{V}_{i\to j}=1\). This constraint suggests that a node can at most transfer the total of its features to other nodes. First, we define the discretized divergence from Equation (2a), that operates on the learned edge weights \(\mathbf{V}\):

\[\mathbf{DIV}_{i}(\mathbf{V}\mathbf{U})=\sum_{j\in\mathcal{N}_{i}}\mathbf{V}_{ j\to i}\odot\mathbf{U}_{j}-\mathbf{U}_{i}\odot\sum_{j\in\mathcal{N}_{i}} \mathbf{V}_{i\to j}=\sum_{j\in\mathcal{N}_{i}}\mathbf{V}_{j\to i}\odot \mathbf{U}_{j}-\mathbf{U}_{i}, \tag{4}\]

where \(\odot\) is the elementwise Hadamard product. Then, the graph advection operator in Algorithm 1 is:

\[\mathbf{U}_{i}^{(l+1/3)}=\mathbf{U}_{i}^{(l)}+h\mathbf{DIV}_{i}(\mathbf{V}^{(l )}\mathbf{U}^{(l)})=\mathbf{U}_{i}^{(l)}+h\left(\sum_{j\in\mathcal{N}_{i}} \mathbf{V}_{j\to i}^{(l)}\odot\mathbf{U}_{j}^{(l)}-\mathbf{U}_{i}^{(l)}\right). \tag{5}\]

Namely, the updated node features are obtained by adding the \(\mathbf{V}_{j\to i}\) weighted inbound node features, while removing the \(\mathbf{V}_{i\to j}\) weighted outbound node features, and \(h\) is a positive step size. The scheme in Equation (5) is the forward Euler discretization. We now show that the proposed graph neural advection operator is _mass conserving_, _stable_ and _consistent_2. By satisfying these properties, our advection operator is adherent to the continuous advection PDE [48].

Footnote 2: See stability definition and proofs in Appendix B.

**Lemma 1**.: _Define the mass of the graph node features \(\mathbf{U}^{(l)}\in\mathbb{R}^{n\times c}\) as the scalar \(\rho^{(l)}=\sum\mathbf{U}^{(l)}\). Then the advection operator in Equation (5) is mass conserving, i.e., \(\rho^{(l+1/3)}=\rho^{(l)}\)._

**Lemma 2**.: _The advection operator in Equation (5) is stable._

To learn a _consistent_ advection operator, that is, an operator that mimics the directional behavior of the advection in Equation (1), we craft an edge weight \(\mathbf{V}\) learning mechanism, shown in Algorithm 2, that yields direction-oriented weights, i.e., we ensure that \(\mathbf{V}_{i\to j}\neq\mathbf{V}_{j\to i}\), unless they are zeroes.

**Input:** Node features \(\mathbf{U}^{(l)}\in\mathbb{R}^{n\times c}\)

**Output:** Edge weights \(\mathbf{V}_{i\to j}^{(l)},\mathbf{V}_{j\to i}^{(l)}\in\mathbb{R}^{c}\)

\(\mathbf{Z}_{ij}^{(l)}=\mathrm{ReLU}(\mathbf{U}_{i}^{(l)}\mathbf{A}_{1}^{(l)} +\mathbf{U}_{j}^{(l)}\mathbf{A}_{2}^{(l)})\mathbf{A}_{3}^{(l)}\quad\mathrm{ and}\quad\mathbf{Z}_{ji}^{(l)}=\mathrm{ReLU}(\mathbf{U}_{j}^{(l)}\mathbf{A}_{1}^{(l)}+ \mathbf{U}_{i}^{(l)}\mathbf{A}_{2}^{(l)})\mathbf{A}_{3}^{(l)}\)

\(\mathbf{V}_{i\to j}^{(l)}=\mathrm{ReLU}(\mathbf{Z}_{ij}^{(l)}-\mathbf{Z}_{ji} ^{(l)})\mathbf{A}_{4}^{(l)}\quad\mathrm{and}\quad\mathbf{V}_{j\to i}^{(l)}= \mathrm{ReLU}(-\mathbf{Z}_{ij}^{(l)}+\mathbf{Z}_{ji}^{(l)})\mathbf{A}_{4}^{(l)}\)

\(\mathbf{V}_{i\to j}^{(l)}\leftarrow\frac{\exp(\mathbf{V}_{i\to j}^{(l)})}{\sum_ {k\in\mathcal{N}_{i}}\exp(\mathbf{V}_{i\to k}^{(l)})}\quad\mathrm{and}\quad \mathbf{V}_{j\to i}^{(l)}\leftarrow\frac{\exp(\mathbf{V}_{j\to i}^{(l)})}{ \sum_{k\in\mathcal{N}_{j}}\exp(\mathbf{V}_{j\to k}^{(l)})}\)

**Algorithm 2** Learning directional edge weights.

Here, \(\boldsymbol{\theta}_{a}^{(l)}=\{\mathbf{A}_{1}^{(l)},\mathbf{A}_{2}^{(l)}, \mathbf{A}_{3}^{(l)},\mathbf{A}_{4}^{(l)}\}\) are learnable fully connected layers, and the \(\exp\) is computed channel-wise. We note that the sign of \(\mathbf{Z}_{ij}-\mathbf{Z}_{ji}\) is opposite than that of \(-\mathbf{Z}_{ij}+\mathbf{Z}_{ji}\) in Algorithm 2. Hence, after the \(\mathrm{ReLU}(\cdot)\) activation, one of the edge weights, either \(\mathbf{V}_{i\to j}\) or \(\mathbf{V}_{j\to i}\) is guaranteed to be equal to zero, and the other will be non-negative. This allows the architecture to create significant asymmetry in the edge weights \(\mathbf{V}\), as also seen in Figure 1.

**Diffusion.** To discretize the diffusion term from Equation (2a), both explicit and implicit time discretizations can be used [4]. An explicit forward Euler discretization yields the following layer:

\[\mathbf{U}^{(l+2/3)}=\mathbf{U}^{(l+1/3)}-h\left(\hat{\mathbf{L}}\mathbf{U}^{(l +1/3)}\mathbf{K}^{(l)}\right). \tag{6}\]

However, an explicit scheme requires using a small step size \(h>0\), as it is marginally stable [4]. We therefore harness an implicit scheme, which guarantees the stability of the diffusion 3, and reads:
\[\mathbf{U}^{(l+2/3)}=\mathrm{mat}\left((\mathbf{I}+h\mathbf{K}^{(l)}\otimes\hat{ \mathbf{L}})^{-1}\mathrm{vec}(\mathbf{U}^{(l+1/3)})\right). \tag{7}\]

Here, \(\otimes\) is the Kronecker product, \(\mathrm{vec}()\) is a flattening operator, and \(\mathrm{mat}()\) reshapes a vector to a matrix. The computation of \(\mathbf{U}^{(l+2/3)}\) requires the solution of a linear system, solved by conjugate gradients4[31, 4]. In our experiments we found 5 iterations to be sufficient.

Footnote 4: We note that the matrix \(\mathbf{I}+h\mathbf{K}_{l}\otimes\hat{\mathbf{L}}\) is positive definite and invertible, because the identity matrix is positive definite, \(h\) is positive, \(\mathbf{K}_{l}\) is non-negative, and the graph Laplacian \(\hat{\mathbf{L}}\) is positive semi-definite.

**Reaction.** Our reaction term is realized using MLPs. Recent works showed that utilizing both additive and multiplicative MLPs yields improved performance [42, 22, 7]. Hence, we define:

\[f(\mathbf{U}^{(l+2/3)},\mathbf{U}^{(0)},\mathbf{X},t;\mathbf{\theta}^{(l)}_{r})= \sigma\left(\mathbf{U}^{(l+2/3)}\mathbf{R}^{(l)}_{1}+\mathrm{tanh}(\mathbf{U }^{(l+2/3)}\mathbf{R}^{(l)}_{2})\odot\mathbf{U}^{(l+2/3)}+\mathbf{U}^{(0)} \mathbf{R}^{(l)}_{3}\right), \tag{8}\]

as our reaction term in Equation (2a). Here, \(\mathbf{\theta}^{(l)}_{r}=\{\mathbf{R}^{(l)}_{1},\mathbf{R}^{(l)}_{2},\mathbf{R} ^{(l)}_{3}\}\) are trainable fully-connected layers, and \(\sigma\) is non-linear activation function (ReLU on our experiments), that can also be coupled with batch-normalization. The reaction term is integrated via forward Euler, as shown in Algorithm 1.

## 4 Experimental Results

We demonstrate our ADR-GNN on two types of tasks on real-world datasets: node classification, and spatio-temporal node forecasting. Architectures and training details are provided in Appendix C, and the complexity of ADR-GNN is discussed in Appendix D. We use a grid search to select hyperparameters, discussed in Appendix E. Datasets details and statistics are reported in Appendix F. Throughout all the experiments, the top three models according to the mean metric are highlighted by **First**, **Second**, **Third**. Overall, we propose the following ADR-GNN architectures:

* ADR-GNN\({}_{S}\). Here we follow a similar approach to typical neural networks, where different weights are learned for each layer. From a dynamical system perspective, this can be interpreted as an unrolled ADR iteration [57]. This architecture is suitable for'static' datasets that do not involve temporal information, such as Cora [59], and is specified in Appendix C.1.
* ADR-GNN\({}_{T}\). A time-dependent ADR-GNN that is suitable for temporal datasets. In addition to the functionality of ADR-GNN\({}_{I}\), it utilizes temporal embedding, discussed in Appendix C.2.

### Node Classification

**Homophilic graphs.** We experiment with Cora [59], Citeseer [76], and Pubmed [60] datasets. We use the 10 splits from [65] with train/validation/test split ratios of \(48\%/32\%/20\%\), and report their average accuracy in Table 1. In Appendix G.1 we also provide the accuracy standard deviation. As a comparison, we consider multiple recent methods, such as GCN [45], GAT [88], Geom-GCN [65], APPNP [46], JKNet [94], MixHop [1], WRGAT[82], GCNII [18], PDE-GCN [27], NSD [10], H2GCN [102], GGCN [95], C&S [40], DMP [96], GREAD [22], LINKX [51], and ACMII [55]. We see that our ADR-GNNs outperforms all methods on the Cora and Pubmed datasets, and achieves close (0.12% accuracy difference) to the best performing PDE-GCN on Citeseer.

**Heterophilic graphs.** While our ADR-GNN offers competitive accuracy on homophilic datasets, as discussed in Section 2, ADR systems are widely used to model non-smooth phenomena and patterns, as often appear in heterophilic datasets by their definition [65]. We therefore utilize 10 heterophilic datasets from various sources. In Table 2 we compare the average accuracy of our ADR-GNN\({}_{\mathrm{S}}\) with recent GNNs on the Squirrel, Film, and Chameleon from [67], as well as the Cornell, Texas and Wisconsin datasets from [65], using the 10-splits from [65] we train/validation/test split ratios of \(48\%/32\%/20\%\). We include more comparisons and the accuracy standard deviation in Appendix G.1. In addition to the previously considered methods, we also compare with FAGCN [9], GraphCON [70], GPR-GNN [21], GRAFF [26], ACMP-GCN [90], and G\({}^{2}\)[71]. We see that ADR-GNNs offers accuracy that is in line with recent state-of-the-art methods. In addition, we evaluate ADR-GNNs on the Twitch-DE, deezer-europe, Penn94, and arXiv-year datasets from [52, 51] to further demonstrate the efficacy of our method, in Appendix G.2.
### Spatio-Temporal Node Forecasting

Classical ARD models are widely utilized to predict and model spatio-temporal phenomena [29; 2]. We therefore now evaluate our temporal ADR-GNN\({}_{\text{T}}\) on several spatio-temporal node forecasting datasets. To this end, we harness the software package PyTorch-Geometric-Temporal [68] that offers a graph machine learning pipeline for spatio-temporal graph tasks. In our experiments, we use the Chickenpox Hungary, PedalMe London, and Wikipedia Math datasets from [68], as well as the traffic speed prediction datasets METR-LA [41] and PEMS-BAY [16].

For the first three datasets, we follow the incremental training mode, mean-squared-error (MSE) loss, and testing procedure from [68]. We report the performance of ADR-GNN\({}_{\text{T}}\) and other models, in terms of MSE, in Table 3. We compare with several recent methods, namely, DCRNN [49], GConv [77], GC-LSTM [17], DyGrAE [85; 84], EGCN [64], A3T-GCN [101], T-GCN [99], MPNN LSTM [63], and AGCRN [5]. Our results in Table 3 show improvement over the considered models, further revealing the significance of neural ADR systems on graphs, offered by our ADR-GNN\({}_{\text{T}}\).

On the METR-LA and PEMS-BAY datsets, we follow the same training and testing procedures, and mean-absolute-error (MAE) loss as in [49]. We report the MAE, root mean squared error (RMSE), and mean absolute percentage error (MAPE). To demonstrate the effectiveness of ADR-GNN\({}_{\text{T}}\) for varying time frame predictions, we report the results on 3, 6, and 12 future frame traffic speed prediction, where each time frame equates to 5 minutes, in Table 4. We compare ADR-GNN\({}_{\text{T}}\) with various methods, from 'classical' (non-neural) approaches such as historical averaging (HA), VAR [54], and SVR [80], to neural methods like FC-LSTM [83], DCRNN [49], Graph WaveNet [93], ASTGCN [34], STSGCN [81], GMAN [100], MTGNN [92], GTS [78], and STEP [79]. We find that our ADR-GNN\({}_{\text{T}}\) offers lower (better) metrics than the considered methods. For instance, on METR-LA, ADR-GNN\({}_{\text{T}}\) reduces the MAE achieved by the recent STEP method from 3.37 to 3.19.

### Ablation Studies

**Synthetic Feature Transportation.** The benefit of diffusion and reaction are known in GNNs (see [30; 14; 22] and references within). However, the significance of neural advection was not studied in GNNs prior to our work, to the best of our knowledge. Therefore, and following the

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline Method & Cora & Cite. & Pub. & & & \\ Homophily & 0.81 & 0.80 & 0.74 & & & \\ \hline GCN & 85.77 & 73.68 & 88.13 & GCN & 23.96 & 26.86 & 28.18 & 52.70 & 52.16 & 48.92 \\ GAT & 86.37 & 74.32 & 87.62 & GAT & 30.03 & 28.45 & 42.93 & 54.32 & 58.38 & 49.41 \\ GCNII\({}^{\dagger}\) & 88.49 & 77.13 & **90.30** & GCNII\({}^{\dagger}\) & 38.47 & 32.87 & 60.61 & 74.86 & 69.46 & 74.12 \\ Geom-GCN\({}^{\dagger}\) & 85.27 & **77.99** & 90.05 & GEOM & 38.32 & 31.63 & 60.90 & 60.81 & 67.57 & 64.12 \\ APPNP & 87.87 & 76.53 & 89.40 & PDE-GCN\({}^{\dagger}\) & - & - & - & - & - \\ JKNet & 85.25 & 75.85 & 88.94 & PDE-GCN\({}^{\dagger}\) & - & - & - & - \\ MixHop & 87.61 & 76.26 & 85.31 & GCN & 23.96 & 26.86 & 85.94 & 85.95 & 89.41 \\ WRGAT & 88.20 & 76.81 & 88.52 & GCN & 23.96 & GCN & 23.96 & 84.86 & 86.86 \\ PDE-GCN\({}^{\dagger}\) & **88.60** & **78.48** & 89.93 & GCN & 23.96 & GCN & 23.96 & 87.65 & 87.65 \\ NSD\({}^{\dagger}\) & 87.14 & 77.14 & 89.49 & GGCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ GGCN & 87.95 & 77.14 & 89.15 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ H2GCN & 87.87 & 77.11 & 89.49 & GGCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ C\&S & **89.05** & 76.22 & 89.74 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ GRAFF\({}^{\dagger}\) & 88.01 & 77.30 & 90.04 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ DMP\({}^{\dagger}\) & 86.52 & 76.87 & 89.27 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ GREAD\({}^{\dagger}\) & 88.57 & 77.60 & **90.23** & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ LINKX & 84.64 & 73.19 & 87.86 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ ACMII\({}^{\dagger}\) & 88.25 & 77.12 & 89.71 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 & GCN & 23.96 \\ \hline ADR-GNN\({}_{\text{S}}\) & **89.43** & **78.36** & **90.55** & GCN & 23.96 & GCN & 23.
discussion of the task of feature transportation in Section 1 and Figure 1, we now compare the behavior of the advection, diffusion, and reaction terms on this task. Although this experiment is conceptually simple, it is evident from Figure 2, that diffusion and reaction terms in GNNs are limited in modeling such a behavior. This result, however, is not surprising. Employing diffusion smooths, rather than directly _transferring_ node features. Similarly, the reaction term can only learn to scale the node features in this experiment. On the contrary, Figure 2 that the advection term, that by definition, transports information, achieves an exact fit. More details about the experiment are given in Appendix G.3.

**The Impact of Advection, Diffusion, and Reaction.** We study the influence of each of the proposed terms in Equation (2a) on real-world datasets, independently and jointly. The results, reported in Table 5 further show the significance of the advection term. For datasets that are homophilic like Cora, we see minor accuracy improvement when incorporating the advection term. This is in line with known findings regarding the benefits of diffusion for homophilic datasets [30; 14]. More importantly, we see that mostly for heterophilic datasets like Chameleon, as well as traffic prediction datasets like PEMS-BAY, utilizing the advection significantly improves the performance of the network.

**The Influence of Number of Layers.** The design of ADR-GNN can alleviate oversmoothing in two ways. First, by learning the diffusion coefficients \(\mathbf{K}\), ADR-GNN controls the amount of smoothing, and can also achieve no smoothing if \(\mathbf{K}\) is zero, depending on the data. Second, note that the advection and reaction terms can increase the frequency of the node features, because they are not limited to smoothing processes. To verify our observation, we evaluate ADR-GNN\({}_{\mathrm{S}}\) on Cora and Citeseer with 2 to 64 layers, to see if its performance degrades as more layers are added, an issue that is associated with oversmoothing. We report the obtained accuracy in Figure 3, where no performance drop is evident. For reference, we also report the results obtained with GCN [45]. Also, we define and report the measured Dirichlet energy in Appendix G.4, that shows ADR-GNN does not oversmooth.

\begin{table}
\begin{tabular}{c c c c} \hline \hline Method & Chickenpox Hungary & PedalMe London & Wikipedia Math \\ \hline DCRNN & 1.124 \(\pm\) 0.015 & 1.463\(\pm\) 0.019 & **0.679 \(\pm\) 0.020** \\ GConvGRU & 1.128 \(\pm\) 0.011 & 1.622 \(\pm\) 0.032 & **0.657 \(\pm\) 0.015** \\ GConvLSTM & 1.121 \(\pm\) 0.014 & **1.442 \(\pm\) 0.028** & 0.777 \(\pm\) 0.021 \\ GC-LSTM & 1.115 \(\pm\) 0.014 & **1.455 \(\pm\) 0.023** & 0.779 \(\pm\) 0.023 \\ DyGrAE & 1.120 \(\pm\) 0.021 & 1.455 \(\pm\) 0.031 & 0.773 \(\pm\) 0.009 \\ EGCN-H & **1.113 \(\pm\) 0.016** & 1.467 \(\pm\) 0.026 & 0.775 \(\pm\) 0.022 \\ EGCN-O & 1.124 \(\pm\) 0.009 & 1.491 \(\pm\) 0.024 & 0.750 \(\pm\) 0.014 \\ A3T-GCN & **1.114 \(\pm\) 0.008** & 1.469 \(\pm\) 0.027 & 0.781 \(\pm\) 0.011 \\ T-GCN & 1.117 \(\pm\) 0.011 & 1.479 \(\pm\) 0.012 & 0.764 \(\pm\) 0.011 \\ MPNN LSTM & 1.116 \(\pm\) 0.023 & 1.485 \(\pm\) 0.028 & 0.795 \(\pm\) 0.010 \\ AGCRN & 1.120 \(\pm\) 0.010 & 1.469 \(\pm\) 0.030 & 0.788 \(\pm\) 0.011 \\ \hline ADR-GNN\({}_{\mathrm{T}}\) & **0.817 \(\pm\) 0.012** & **0.598 \(\pm\) 0.050** & **0.571 \(\pm\) 0.014** \\ \hline \hline \end{tabular}
\end{table}
Table 3: The predictive performance of spatio-temporal neural networks evaluated by average MSE of 10 experimental repetitions and standard deviations, calculated on 10% forecasting horizons.

Figure 2: Source and target node features, and their fit using advection, diffusion, and reaction.
## 5 Limitations

While our ADR-GNN allows modeling data that involves directed transportation of features, it comes with an additional cost, as it requires computing the directed edge weights to mimic the behavior of the advection PDE. This adds complexity to each layer in the form of MLPs, described in Algorithm 2. Furthermore, as we show in our experiments, incorporating the advection into GNNs yields minor improvement when the data are governed by diffusion, e.g., for the Pubmed dataset.
Summary and Discussion

In this paper, we present a novel GNN architecture that is based on the Advection-Diffusion-Reaction PDE, called ADR-GNN. We develop a graph neural advection operator that mimics the continuous advection operator, and compose it with learnable diffusion and reaction terms.

We discuss and analyze the properties of ADR-GNN and its flexibility in modeling various phenomena. In particular, we show that the main advantage of the graph advection operator is its ability to transport information over the graph edges through the layers - a behavior that is hard to model using the diffusion and reaction terms that have been used in the literature. To demonstrate the effectiveness of ADR-GNN we experiment with total of 18 real-world datasets, from homophilic and heterophilic node classification to spatio-temporal node forecasting datasets.

While the gains observed on homophilic datasets are relatively modest, the performance improvements demonstrated on heterophilic datasets are significant, offering 5% accuracy increase in some cases. Moreover, when applied to spatio-temporal node forecasting datasets, our ADR-GNN exhibits notable enhancements in the evaluated metrics compared to other methods. This progress can be attributed to the inherent suitability of ADR-GNN for tasks involving directional transportation of features, making it an intuitive choice for modeling such scenarios.

## References

* Abu-El-Haija et al. [2019] Sami Abu-El-Haija, Bryan Perozzi, Amol Kapoor, Nazanin Alipourfard, Kristina Lerman, Hrayr Harutyunyan, Greg Ver Steeg, and Aram Galstyan. Mixhop: Higher-order graph convolutional architectures via sparsified neighborhood mixing. In _international conference on machine learning_, pages 21-29. PMLR, 2019.
* Adam and Sibert [2004] M Shiham Adam and John R Sibert. Use of neural networks with advection-diffusion-reaction models to estimate large-scale movements of skipjack tuna from tagging data. Technical report, Marine Research Centre, Ministry of Fisheries, Agriculture and Marine Resources, 2004.
* Alon and Yahav [2021] Uri Alon and Eran Yahav. On the bottleneck of graph neural networks and its practical implications. In _International Conference on Learning Representations_, 2021.
* Ascher [2008] Uri M Ascher. _Numerical methods for evolutionary differential equations_. SIAM, 2008.
* Bai et al. [2020] Lei Bai, Lina Yao, Can Li, Xianzhi Wang, and Can Wang. Adaptive Graph Convolutional Recurrent Network for Traffic Forecasting. _Advances in Neural Information Processing Systems_, 33, 2020.
* de Avila Belbute-Peres et al. [2020] Filipe de Avila Belbute-Peres, Thomas D. Economon, and J. Zico Kolter. Combining Differentiable PDE Solvers and Graph Neural Networks for Fluid Flow Prediction. In _International Conference on Machine Learning (ICML)_, 2020.
* Ben-Shaul et al. [2023] Ido Ben-Shaul, Tomer Galanti, and Shai Dekel. Exploring the approximation capabilities of multiplicative neural networks for smooth functions. _arXiv preprint arXiv:2301.04605_, 2023.
* Betts [2001] J.T. Betts. _Practical Methods for Optimal Control using Nonlinear Programming_. Advances in Design and Control. SIAM, Philadelphia, 2001.
* Bo et al. [2021] Deyu Bo, Xiao Wang, Chuan Shi, and Huawei Shen. Beyond low-frequency information in graph convolutional networks. In _AAAI_. AAAI Press, 2021.
* Bodnar et al. [2022] Cristian Bodnar, Francesco Di Giovanni, Benjamin Paul Chamberlain, Pietro Lio, and Michael M Bronstein. Neural sheaf diffusion: A topological perspective on heterophily and oversmoothing in gnns. _arXiv preprint arXiv:2202.04579_, 2022.
* Borovitskiy et al. [2021] Viacheslav Borovitskiy, Iskander Azangulov, Alexander Terenin, Peter Mostowsky, Marc Peter Deisenroth, and Nicolas Durrande. Matern gaussian processes on graphs. In _International Conference on Artificial Intelligence and Statistics_. PMLR, 2021.
* Brandstetter et al. [2022] Johannes Brandstetter, Daniel E. Worrall, and Max Welling. Message passing neural PDE solvers. In _International Conference on Learning Representations_, 2022.
* Cai and Wang [2020] Chen Cai and Yusu Wang. A note on over-smoothing for graph neural networks. _arXiv preprint arXiv:2006.13318_, 2020.
* [14] Benjamin Paul Chamberlain, James Rowbottom, Maria Gorinova, Stefan Webb, Emanuele Rossi, and Michael M Bronstein. GRAND: Graph neural diffusion. In _International Conference on Machine Learning (ICML)_, pages 1407-1418. PMLR, 2021.
* [15] Airlie Chapman and Airlie Chapman. Advection on graphs. _Semi-Autonomous Networks: Effective Control of Networked Systems through Protocols, Design, and Modeling_, pages 3-16, 2015.
* [16] Chao Chen, Karl Petty, Alexander Skabardonis, Pravin Varaiya, and Zhanfeng Jia. Freeway performance measurement system: mining loop detector data. _Transportation Research Record_, 1748(1):96-102, 2001.
* [17] Jinyin Chen, Xuanheng Xu, Yangyang Wu, and Haibin Zheng. GC-LSTM: Graph Convolution Embedded LSTM for Dynamic Link Prediction. _arXiv preprint arXiv:1812.04206_, 2018.
* [18] Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li. Simple and deep graph convolutional networks. In Hal Daume III and Aarti Singh, editors, _Proceedings of the 37th International Conference on Machine Learning_, volume 119 of _Proceedings of Machine Learning Research_, pages 1725-1735. PMLR, 13-18 Jul 2020.
* [19] Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K. Duvenaud. Neural ordinary differential equations. _CoRR_, abs/1806.07366, 2018.
* [20] Weikai Chen, Xiaoguang Han, Guanbin Li, Chao Chen, Jun Xing, Yajie Zhao, and Hao Li. Deep rbfnet: Point cloud feature learning using radial basis functions. _arXiv preprint arXiv:1812.04302_, 2018.
* [21] Eli Chien, Jianhao Peng, Pan Li, and Olgica Milenkovic. Adaptive universal generalized pagerank graph neural network. In _International Conference on Learning Representations_, 2021.
* [22] Jeongwhan Choi, Seoyoung Hong, Noseong Park, and Sung-Bae Cho. Gread: Graph neural reaction-diffusion equations. _arXiv preprint arXiv:2211.14208_, 2022.
* [23] Keith H Coats. A note on impes and some impes-based simulation models. _SPE Journal_, 5(03):245-251, 2000.
* [24] Chris Cosner. Reaction-diffusion-advection models for the effects and evolution of dispersal. _Discrete & Continuous Dynamical Systems_, 34(5):1701, 2014.
* [25] M Deepa Maheshvare, Soumyendu Raha, and Debnath Pal. A graph-based framework for multiscale modeling of physiological transport. _Frontiers in Network Physiology_, 1:18, 2022.
* [26] Francesco Di Giovanni, James Rowbottom, Benjamin P Chamberlain, Thomas Markovich, and Michael M Bronstein. Graph neural networks as gradient flows. _arXiv preprint arXiv:2206.10991_, 2022.
* [27] Moshe Eliasof, Eldad Haber, and Eran Treister. PDE-GCN: Novel architectures for graph neural networks motivated by partial differential equations. _Advances in Neural Information Processing Systems_, 34:3836-3849, 2021.
* [28] L. C. Evans. _Partial Differential Equations_. American Mathematical Society, San Francisco, 1998.
* [29] Bernold Fiedler and Arnd Scheel. Spatio-temporal dynamics of reaction-diffusion patterns. _Trends in nonlinear analysis_, pages 23-152, 2003.
* [30] Johannes Gasteiger, Stefan Weissenberger, and Stephan Gunnemann. Diffusion improves graph learning. In _Conference on Neural Information Processing Systems (NeurIPS)_, 2019.
* [31] G.H. Golub and C.F. van Loan. _Matrix Computations_. Johns Hopkins University Press, 1988.
* [32] Alessio Gravina, Davide Bacciu, and Claudio Gallicchio. Anti-symmetric dgn: a stable architecture for deep graph networks. _arXiv preprint arXiv:2210.09789_, 2022.
* [33] Fangda Gu, Heng Chang, Wenwu Zhu, Somayeh Sojoudi, and Laurent El Ghaoui. Implicit graph neural networks. _Advances in Neural Information Processing Systems_, 33:11984-11995, 2020.
* [34] Shengnan Guo, Youfang Lin, Ning Feng, Chao Song, and Huaiyu Wan. Attention based spatial-temporal graph convolutional networks for traffic flow forecasting. In _Proceedings of the AAAI conference on artificial intelligence_, volume 33, pages 922-929, 2019.
* [35] E. Haber. _Computational Methods in Geophysical Electromagnetics_. SIAM, Philadelphia, 2014.
* [36] E. Haber and L. Ruthotto. Stable architectures for deep neural networks. _Inverse Problems_, 34(1), 2017.
* [37] Eldad Haber, Keegan Lensink, Eran Treister, and Lars Ruthotto. IMEXnet a forward stable deep neural network. In _International Conference on Machine Learning_, pages 2525-2534. PMLR, 2019.
* [38] M. Hochbruck, C. Luibich, and H. Selhofer. Exponential integrators for large systems of differential equations. _SIAN. J. Sci. Comp_, 19, 1998. n5.
* [39] Radim Hosek and Jonas Volek. Discrete advection-diffusion equations on graphs: Maximum principle and finite volumes. _Applied Mathematics and Computation_, 361:630-644, 2019.
* [40] Qian Huang, Horace He, Abhay Singh, Ser-Nam Lim, and Austin R Benson. Combining label propagation and simple models out-performs graph neural networks. _arXiv preprint arXiv:2010.13993_, 2020.
* [41] Hosagrahar V Jagadish, Johannes Gehrke, Alexandros Labrinidis, Yannis Papakonstantinou, Jignesh M Patel, Raghu Ramakrishnan, and Cyrus Shahabi. Big data and its technical challenges. _Communications of the ACM_, 57(7):86-94, 2014.
* [42] Siddhant M. Jayakumar, Wojciech M. Czarnecki, Jacob Menick, Jonathan Schwarz, Jack Rae, Simon Osindero, Yee Whye Teh, Tim Harley, and Razvan Pascanu. Multiplicative interactions and where to find them. In _International Conference on Learning Representations_, 2020.
* [43] Samet Kadioglu, Dana Knoll, Mark Sussman, and Richard Martineau. A second order jfnk-based imex method for single and multi-phase flows. In _Computational Fluid Dynamics 2010_, pages 549-554. Springer, 2011.
* [44] Alexei M Khokhlov. Propagation of turbulent flames in supernovae. _The Astrophysical Journal_, 449:695, 1995.
* [45] Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. _International Conference on Learning Representations (ICLR)_, 2017.
* [46] Johannes Klicpera, Aleksandar Bojchevski, and Stephan Gunnemann. Combining neural networks with personalized pagerank for classification on graphs. In _International Conference on Learning Representations_, 2019.
* [47] Toshiyuki Koto. Imex runge-kutta schemes for reaction-diffusion equations. _Journal of Computational and Applied Mathematics_, 215(1):182-195, 2008.
* [48] R.J. LeVeque. _Numerical Methods for Conservation Laws_. Birkhauser, 1990.
* [49] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion Convolutional Recurrent Neural Network: Data-Driven Traffic Forecasting. In _International Conference on Learning Representations_, 2018.
* [50] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. _Advances in Neural Information Processing Systems_, 33:6755-6766, 2020.
* [51] Derek Lim, Felix Matthew Hohne, Xiuyu Li, Sijia Linda Huang, Vaishnavi Gupta, Omkar Prasad Bhalerao, and Ser-Nam Lim. Large scale learning on non-homophilous graphs: New benchmarks and strong simple methods. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_, 2021.
* [52] Derek Lim, Xiuyu Li, Felix Hohne, and Ser-Nam Lim. New benchmarks for learning on non-homophilous graphs. _Workshop on Graph Learning Benchmarks, WWW_, 2021.
* [53] Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. PDE-net: Learning PDEs from data, 2018.
* [54] Zheng Lu, Chen Zhou, Jing Wu, Hao Jiang, and Songyue Cui. Integrating granger causality and vector auto-regression for traffic prediction of large-scale wlans. _KSII Transactions on Internet and Information Systems (TIIS)_, 10(1):136-151, 2016.
* [55] Sitao Luan, Chenqing Hua, Qincheng Lu, Jiaqi Zhu, Mingde Zhao, Shuyuan Zhang, Xiao-Wen Chang, and Doina Precup. Revisiting heterophily for graph neural networks. _Conference on Neural Information Processing Systems_, 2022.
* [56] Sitao Luan, Mingde Zhao, Chenqing Hua, Xiao-Wen Chang, and Doina Precup. Complete the missing half: Augmenting aggregation filtering with diversification for graph convolutional networks. _arXiv preprint arXiv:2008.08844_, 2020.
* [57] Morteza Mardani, Qingyun Sun, David Donoho, Vardan Papyan, Hatef Monajemi, Shreyas Vasanawala, and John Pauly. Neural proximal gradient descent for compressive imaging. _Advances in Neural Information Processing Systems_, 31, 2018.
* [58] Robert MM Mattheij, Sjoerd W Rienstra, and JHM Ten Thije Boonkkamp. _Partial differential equations: modeling, analysis, computation_. SIAM, 2005.
* [59] Andrew Kachites McCallum, Kamal Nigam, Jason Rennie, and Kristie Seymore. Automating the construction of internet portals with machine learning. _Information Retrieval_, 3(2):127-163, 2000.
* [60] Galileo Namata, Ben London, Lise Getoor, Bert Huang, and U Edu. Query-driven active surveying for collective classification. In _10th International Workshop on Mining and Learning with Graphs_, volume 8, page 1, 2012.
* [61] Hoang Nt and Takanori Maehara. Revisiting graph neural networks: All we have is low-pass filters. _arXiv preprint arXiv:1905.09550_, 2019.
* [62] Kenta Oono and Taiji Suzuki. Graph neural networks exponentially lose expressive power for node classification. In _International Conference on Learning Representations_, 2020.
* [63] George Panagopoulos, Giannis Nikolentzos, and Michalis Vazirgiannis. Transfer Graph Neural Networks for Pandemic Forecasting. In _Proceedings of the 35th AAAI Conference on Artificial Intelligence_, 2021.
* [64] Aldo Pareja, Giacomo Domeniconi, Jie Chen, Tengfei Ma, Toyotaro Suzumura, Hiroki Kanezashi, Tim Kaler, Tao B Schardl, and Charles E Leiserson. EvolveGCN: Evolving Graph Convolutional Networks for Dynamic Graphs. In _AAAI_, pages 5363-5370, 2020.
* [65] Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang. Geom-gcn: Geometric graph convolutional networks. In _International Conference on Learning Representations_, 2020.
* [66] Maziar Raissi. Deep hidden physics models: Deep learning of nonlinear partial differential equations. _The Journal of Machine Learning Research_, 19(1):932-955, 2018.
* [67] Benedek Rozemberczki, Carl Allen, and Rik Sarkar. Multi-Scale Attributed Node Embedding. _Journal of Complex Networks_, 9(2), 2021.
* [68] Benedek Rozemberczki, Paul Scherer, Yixuan He, George Panagopoulos, Alexander Riedel, Maria Astefanoaei, Oliver Kiss, Ferenc Beres, Guzman Lopez, Nicolas Collignon, et al. Pytorch geometric temporal: Spatiotemporal signal processing with neural machine learning models. In _Proceedings of the 30th ACM International Conference on Information & Knowledge Management_, pages 4564-4573, 2021.
* [69] T Konstantin Rusch, Michael M Bronstein, and Siddhartha Mishra. A survey on oversmoothing in graph neural networks. _arXiv preprint arXiv:2303.10993_, 2023.
* [70] T Konstantin Rusch, Ben Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael Bronstein. Graph-coupled oscillator networks. In _International Conference on Machine Learning_, pages 18888-18909. PMLR, 2022.
* [71] T Konstantin Rusch, Benjamin P Chamberlain, Michael W Mahoney, Michael M Bronstein, and Siddhartha Mishra. Gradient gating for deep multi-rate learning on graphs. _arXiv preprint arXiv:2210.00513_, 2022.
* [72] Lars Ruthotto and Eldad Haber. Deep neural networks motivated by partial differential equations. _Journal of Mathematical Imaging and Vision_, 62:352-364, 2020.
* [73] MH Saadat, B Gjorgiev, L Das, and G Sansavini. Neural tangent kernel analysis of pinn for advection-diffusion equation. _arXiv preprint arXiv:2211.11716_, 2022.
* [74] Alvaro Sanchez-Gonzalez, Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian graph networks with ode integrators. _arXiv preprint arXiv:1909.12790_, 2019.
* [75] AW Seed. A dynamic and spatial scaling approach to advection forecasting. _Journal of Applied Meteorology and Climatology_, 42(3):381-388, 2003.
* [76] Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Eliassi-Rad. Collective classification in network data. _AI magazine_, 29(3):93-93, 2008.
* [77] Youngjoo Seo, Michael Defferrard, Pierre Vandergheynst, and Xavier Bresson. Structured Sequence Modeling with Graph Convolutional Recurrent Networks. In _International Conference on Neural Information Processing_, pages 362-373. Springer, 2018.
* [78] Chao Shang, Jie Chen, and Jinbo Bi. Discrete graph structure learning for forecasting multiple time series. In _International Conference on Learning Representations_, 2021.
* [79] Zezhi Shao, Zhao Zhang, Fei Wang, and Yongjun Xu. Pre-training enhanced spatial-temporal graph neural network for multivariate time series forecasting. In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_, pages 1567-1577, 2022.
* [80] Alex J Smola and Bernhard Scholkopf. A tutorial on support vector regression. _Statistics and computing_, 14:199-222, 2004.
* [81] Chao Song, Youfang Lin, Shengnan Guo, and Huaiyu Wan. Spatial-temporal synchronous graph convolutional networks: A new framework for spatial-temporal network data forecasting. In _Proceedings of the AAAI conference on artificial intelligence_, volume 34, pages 914-921, 2020.
* [82] Susheel Suresh, Vinith Budde, Jennifer Neville, Pan Li, and Jianzhu Ma. Breaking the limit of graph neural networks by improving the assortativity of graphs with local mixing patterns. _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, 2021.
* [83] Ilya Sutskever, Oriol Vinyals, and Quoc V Le. Sequence to sequence learning with neural networks. _Advances in neural information processing systems_, 27, 2014.
* [84] Aynaz Taheri and Tanya Berger-Wolf. Predictive Temporal Embedding of Dynamic Graphs. In _Proceedings of the 2019 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining_, pages 57-64, 2019.
* [85] Aynaz Taheri, Kevin Gimpel, and Tanya Berger-Wolf. Learning to represent the evolution of dynamic graphs with recurrent models. In _Companion Proceedings of The 2019 World Wide Web Conference_, WWW '19, page 301-307, 2019.
* [86] Matthew Thorpe, Tan Minh Nguyen, Hedi Xia, Thomas Strohmer, Andrea Bertozzi, Stanley Osher, and Bao Wang. GRAND++: Graph neural diffusion with a source term. In _International Conference on Learning Representations_, 2022.
* [87] Lafras Uys. _Coupling kinetic models and advection-diffusion equations to model vascular transport in plants, applied to sucrose accumulation in sugarcane_. PhD thesis, Stellenbosch: University of Stellenbosch, 2009.
* [88] Petar Velickovic, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph Attention Networks. _International Conference on Learning Representations_, 2018.
* [89] Luminita A Vese and Tony F Chan. A multiphase level set framework for image segmentation using the mumford and shah model. _International journal of computer vision_, 50(3):271-293, 2002.
* [90] Yuelin Wang, Kai Yi, Xinliang Liu, Yu Guang Wang, and Shi Jin. Acmp: Allen-cahn message passing for graph neural networks with particle phase transition. _arXiv preprint arXiv:2206.05437_, 2022.
* [91] E Weinan. A Proposal on Machine Learning via Dynamical Systems. _Communications in Mathematics and Statistics_, 5(1):1-11, March 2017.
* [92] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang. Connecting the dots: Multivariate time series forecasting with graph neural networks. In _Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining_, pages 753-763, 2020.
* [93] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, and Chengqi Zhang. Graph wavenet for deep spatial-temporal graph modeling. In _Proceedings of the 28th International Joint Conference on Artificial Intelligence_, IJCAI'19, page 1907-1913. AAAI Press, 2019.
* [94] Keyulu Xu, Chengtao Li, Yonglong Tian, Tomohiro Sonobe, Ken-ichi Kawarabayashi, and Stefanie Jegelka. Representation learning on graphs with jumping knowledge networks. In Jennifer Dy and Andreas Krause, editors, _Proceedings of the 35th International Conference on Machine Learning_, volume 80 of _Proceedings of Machine Learning Research_, pages 5453-5462. PMLR, 10-15 Jul 2018.
* [95] Yujun Yan, Milad Hashemi, Kevin Swersky, Yaoqing Yang, and Danai Koutra. Two sides of the same coin: Heterophily and oversmoothing in graph convolutional neural networks. _arXiv preprint arXiv:2102.06462_, 2021.
* [96] Liang Yang, Mengzhe Li, Liyang Liu, Chuan Wang, Xiaochun Cao, Yuanfang Guo, et al. Diverse message passing for attribute with heterophily. _Advances in Neural Information Processing Systems_, 34:4751-4763, 2021.
* [97] Kuangen Zhang, Ming Hao, Jing Wang, Clarence W de Silva, and Chenglong Fu. Linked dynamic graph cnn: Learning on point cloud via linking hierarchical features. _arXiv preprint arXiv:1904.10014_, 2019.
* [98] Tianjun Zhang, Zhewei Yao, Amir Gholami, Joseph E Gonzalez, Kurt Keutzer, Michael W Mahoney, and George Biros. ANODEV2: A coupled neural ode framework. _Advances in Neural Information Processing Systems_, 32, 2019.
* [99] Ling Zhao, Yujiao Song, Chao Zhang, Yu Liu, Pu Wang, Tao Lin, Min Deng, and Haifeng Li. T-GCN: A Temporal Graph Convolutional Network for Traffic Prediction. _IEEE Transactions on Intelligent Transportation Systems_, 21(9):3848-3858, 2019.
* [100] Chuanpan Zheng, Xiaoliang Fan, Cheng Wang, and Jianzhong Qi. Gman: A graph multi-attention network for traffic prediction. In _Proceedings of the AAAI conference on artificial intelligence_, volume 34, pages 1234-1241, 2020.
* [101] Jiawei Zhu, Yujiao Song, Ling Zhao, and Haifeng Li. A3T-GCN: Attention Temporal Graph Convolutional Network for Traffic Forecasting. _arXiv preprint arXiv:2006.11583_, 2020.
* [102] Jiong Zhu, Yujun Yan, Lingxiao Zhao, Mark Heimann, Leman Akoglu, and Danai Koutra. Beyond homophily in graph neural networks: Current limitations and effective designs. _Advances in Neural Information Processing Systems_, 33:7793-7804, 2020.
* [103] Juntang Zhuang, Nicha Dvornek, Xiaoxiao Li, and James S Duncan. Ordinary differential equations on graph networks. 2020.
The Behavior of the Operator Splitting Approach

The advantage of Operator Splitting (OS) is that it allows the individual treatment of each component of the ODE separately, thus obtaining the appropriate qualitative behavior. This is especially beneficial in the context of advection, where it is difficult to obtain stability, that is, to satisfy the CFL condition [48], as well as to obtain mass conservation. Furthermore, OS is beneficial for the diffusion component, where implicit methods guarantee stability [4] compared to unstable explicit discretizations.

We now discuss the behavior of the OS approach for integrating the ADR ODE in Equation (2). The theory behind OS can be analyzed in the linear case, and in the case of non-linear equations, linearization is typically assumed [58]. We now consider a linear ODE of the form

\[\frac{d\mathbf{U}(t)}{dt}=\mathbf{A}\mathbf{U}+\mathbf{D}\mathbf{U}+\mathbf{R }\mathbf{U}, \tag{9}\]

where \(\mathbf{A}\,\mathbf{D}\,\mathbf{R}\) denote the advection, diffusion, and reaction operators, respectively. \(\mathbf{U}(t)\) denotes the node features at time \(t\).

Suppose that we are interested in computing \(\mathbf{U}(t+\delta t)\). For this constant ODE system, the _analytic_, exact solution is given by [28]:

\[\mathbf{U}(t+\delta t)=\exp\left(\delta t(\mathbf{A}+\mathbf{D}+\mathbf{R}) \right)\mathbf{U}(t). \tag{10}\]

The following Lemma is easily proven using Taylor series of the matrix exponential function [4]:

**Lemma 3**.: _Let \(\mathbf{Q}=\exp(\delta t(\mathbf{A}+\mathbf{D}+\mathbf{R}))\) where \(\mathbf{A}\), \(\mathbf{D}\), and \(\mathbf{R}\) are matrices that do not share their eigenvectors. Then the discrepancy between the exact and OS solution operator reads_

\[\mathbf{Q}-\exp(\delta t\mathbf{R})\exp(\delta t\mathbf{D})\exp(\delta t \mathbf{A})=\mathcal{O}(\delta t^{2})\]

**Remark 1**.: _If the eigenvectors of \(\mathbf{A},\mathbf{D},\mathbf{R}\) from Lemma 3 are shared, then the matrix exponents commute, and the discrepancy is zero._

Following Lemma 3, it holds that the solution of the ADR ODE can be expressed as a sequential process of three separate problems, with an error of \(\mathcal{O}(\delta t^{2})\) compared to the exact solution, as follows:

\[\exp\left(\delta t(\mathbf{A}+\mathbf{D}+\mathbf{R})\right) \mathbf{U}(t) = \exp\left(\delta t\mathbf{R}\right)\exp\left(\delta t\mathbf{D} \right)\exp\left(\delta t\mathbf{A}\right)\mathbf{U}(t)+\mathcal{O}(\delta t ^{2})\] \[= \underbrace{\exp(\delta t\mathbf{R})\underbrace{\exp\left(\delta t \mathbf{D}\right)\underbrace{\exp\left(\delta t\mathbf{A}\right)\mathbf{U}(t)} _{\text{Advection}}}_{\text{Advection}-\text{Diffusion}-\text{Reaction}}}_{ \text{Advection}-\text{Diffusion}-\text{Reaction}}+\mathcal{O}(\delta t^{2})\]

Note that \(\mathbf{U}^{(l+1/3)}=\exp(\delta t\mathbf{A})\mathbf{U}^{(l)}\) is the exact solution [38] of the reaction system \(\frac{\mathbf{U}(t)}{dt}=\mathbf{A}\mathbf{U}\). Similarly, \(\mathbf{U}^{(l+2/3)}=\exp(\delta t\mathbf{D})\mathbf{U}^{(l+1/3)}\) is the exact solution of the diffusion system \(\frac{\mathbf{U}(t)}{dt}=\mathbf{D}\mathbf{U}\), and \(\mathbf{U}^{(l+1)}=\exp(\delta t\mathbf{R})\mathbf{U}^{(l+2/3)}\) is the exact solution of the reaction system \(\frac{\mathbf{U}(t)}{dt}=\mathbf{R}\mathbf{U}\). Thus, operator splitting can be viewed as taking a step of advection, followed by a step of diffusion that is finally followed by a reaction step. Note that the error above is of the same magnitude as each step of the forward Euler integration scheme, often used in neural networks, e.g., in ResNet (see [A] in additional references below).

## Appendix B The Properties of the Graph Discretized Advection Operator

In this section, we prove Lemma 1 and Lemma 2 from the main text. For convenience, we repeat the Lemmas, followed by their proofs. We start by noting the following remark:

**Remark 2**.: _The advection operator in Equation (5) does not mix the node feature \(\mathbf{U}\) channels in our ADR-GNN._

The importance of this remark is that it allows us to analyze the properties of the advection operator per-channel, or, alternatively, assuming a single channel, \(c=1\), which we assume in our following proofs.
**Lemma 1**.: _Define the mass of the graph node features \(\mathbf{U}^{(l)}\in\mathbb{R}^{n\times c}\) as the scalar \(\rho^{(l)}=\sum\mathbf{U}^{(l)}\). Then the advection operator in Equation (5) is mass conserving, i.e., \(\rho^{(l+1/3)}=\rho^{(l)}\)._

Proof.: Without loss of generality, and following Remark 2, let us assume a single channel, and consider the mass of the node features \(\rho=\sum\mathbf{U}\), before and after applying an advection layer as described in Equation (5). The input node features has a total mass of \(\rho^{(l)}=\sum\mathbf{U}^{(l)}\). The total mass of the output of an advection layer reads:

\[\rho^{(l+1/3)}=\sum_{i}\mathbf{U}^{(l+1/3)}_{i} = \sum_{i}\left(\mathbf{U}^{(l)}_{i}+h\sum_{j\in\mathcal{N}_{i}} \mathbf{V}^{(l)}_{j\to i}\mathbf{U}^{(l)}_{j}-h\mathbf{U}^{(l)}_{i}\right) \tag{11a}\] \[= \sum_{i}\left(\mathbf{U}^{(l)}_{i}+h\sum_{j}\mathbf{V}^{(l)}_{j \to i}\mathbf{U}^{(l)}_{j}-h\mathbf{U}^{(l)}_{i}\right)\] (11b) \[= \sum_{i}\mathbf{U}^{(l)}_{i}+h\left(\sum_{i}\sum_{j}\mathbf{U}^{ (l)}_{j}\left(\mathbf{V}^{(l)}_{j\to i}\right)-\sum_{i}\mathbf{U}^{(l)}_{i}\right)\] (11c) \[= \sum_{i}\mathbf{U}^{(l)}_{i}+h\left(\sum_{j}\mathbf{U}^{(l)}_{j} \left(\sum_{i}\mathbf{V}^{(l)}_{j\to i}\right)-\sum_{i}\mathbf{U}^{(l)}_{i}\right)\] (11d) \[= \sum_{i}\mathbf{U}^{(l)}_{i}+h\left(\sum_{j}\mathbf{U}^{(l)}_{j} -\sum_{i}\mathbf{U}^{(l)}_{i}\right)=\sum_{i}\mathbf{U}^{(l)}_{i}=\rho^{(l)}. \tag{11e}\]

The transition between Equations (11a) and (11b) is valid because for \(j\notin\mathcal{N}_{i}\), the edge weight is zero, i.e., \(\mathbf{V}^{(l)}_{j\to i}=0\), and therefore the summation does not change. Also, the transition between Equations (11d) and (11e) holds because of the constraint on the sum of outbound edge weights to be equal to \(1\), i.e. \(\sum_{i\in\mathcal{N}_{j}}\mathbf{V}^{(l)}_{j\to i}=1\). Therefore, because \(\rho^{(l+1/3)}=\rho^{(l)}\), our graph discretized advection operator is mass preserving. 

**Definition 1**.: _A neural operator \(F\) that considers node features \(\mathbf{U}\) is (Lyapunov) stable if for every \(\epsilon>0\) there exists \(\delta>0\), such that every pair of inputs \(\mathbf{U}\,,\tilde{\mathbf{U}}\) that satisfy \(\|\mathbf{U}-\tilde{\mathbf{U}}\|\leq\delta\), then \(\|F(\mathbf{U})-F(\tilde{\mathbf{U}})\|\leq\epsilon\)._

**Lemma 2**.: _The advection operator in Equation (5) is stable._

Proof.: Without loss of generality, and following Remark 2, let us consider a single channel, and let \(\mathbf{V}^{(l)}\) be a sparse matrix such that \(\mathbf{V}^{(l)}_{ij}=\mathbf{V}^{(l)}_{i\to j}\). To show stability, we first observe that the matrix form of the scalar formulation of the advection _layer_ in Equation (5) is given by:

\[\mathbf{U}^{(l+1/3)}=\mathbf{I}\mathbf{U}^{(l)}+h\mathbf{V}^{(l)}\mathbf{U}^ {(l)}-h\mathbf{U}^{(l)}=\underbrace{\left((1-h)\mathbf{I}+h\mathbf{V}^{(l)} \right)}_{\mathbf{A}^{(l)}}\mathbf{U}^{(l)}. \tag{12}\]

Because of the demand that \(\mathbf{V}^{(l)}\) is normalized (i.e., \(\sum_{j\in\mathcal{N}_{i}}\mathbf{V}^{(l)}_{i\to j}=1\)), and the advection weights satisfy \(0\leq\mathbf{V}^{(l)}_{i\to j}\leq 1\), the advection operator \(\mathbf{A}^{(l)}\) is a column stochastic non-negative matrix. By the Perron-Frobenius theorem, such matrices are known to have a spectral radius bounded by 1 (see [B] in additional appendix references), and hence the advection operator is stable. 

## Appendix C Architectures and Training Details

As discussed in the main paper, we propose two architectures, depending on the type of dataset - static (e.g., Cora), or spatio-temporal (e.g., PEMS-BAY). In the following subsections, we elaborate on these architectures.
### Node Classification: ADR-GNNs

We now elaborate on the'static' architecture ADR-GNNs used in our node classification experiments. The overall architecture is similar to standard GNN architectures for node classification, such as GCN [45] and GCNII [18]. It is composed of an initial embedding layer (that corresponds to Equation (2b) in the main paper), \(L\) graph neural ADR layers, and a classifier, as described in Equation (3) in the main paper. The complete flow of ADR-GNNs is described in Algorithm 3. To train ADR-GNNs on node classification datasets, we minimize the cross-entropy loss between the ground-truth node labels \(\mathbf{Y}\) and the predicted node labels \(\tilde{\mathbf{Y}}\), as is standard in GNNs, and similar to [45].

```
Input: Node features \(\mathbf{X}\in\mathbb{R}^{n\times c_{in}}\) Output: Predicted node labels \(\tilde{\mathbf{Y}}\in\mathbb{R}^{n\times c_{out}}\)
1:procedureADR-GNNs
2:\(\mathbf{X}\leftarrow\mathrm{Dropout}(\mathbf{X},p)\)
3:\(\mathbf{U}^{(0)}=g_{in}(\mathbf{X})\)
4:for\(l=0\dots L-1\)do
5:\(\mathbf{U}^{(l)}\leftarrow\mathrm{Dropout}(\mathbf{U}^{(l)},p)\)
6: Advection: \(\mathbf{U}^{(l+1/3)}=\mathbf{U}^{(l)}+h\mathbf{D}\mathbf{I}\mathbf{V}(\mathbf{ V}(\mathbf{U}^{(l)};\boldsymbol{\theta}_{a}^{(l)})\mathbf{U}^{(l)})\)
7: Diffusion: \(\mathbf{U}^{(l+2/3)}=\mathrm{mat}\left((\mathbf{I}+h\mathbf{K}(\boldsymbol{ \theta}_{d}^{(l)})\otimes\hat{\mathbf{L}})^{-1}\mathrm{vec}(\mathbf{U}^{(l+1 /3)})\right)\)
8: Reaction: \(\mathbf{U}^{(l+1)}=\mathbf{U}^{(l+2/3)}+hf(\mathbf{U}^{(l+2/3)},\mathbf{U}^{(0) };\boldsymbol{\theta}_{r}^{(l)})\)
9:endfor
10:\(\mathbf{U}^{(L)}\leftarrow\mathrm{Dropout}(\mathbf{U}^{(L)},p)\)
11:\(\tilde{\mathbf{Y}}=g_{out}(\mathbf{U}^{(L)})\)
12: Return \(\tilde{\mathbf{Y}}\)
13:endprocedure
```

**Algorithm 3** ADR-GNNs Architecture Flow

### Spatio-Temporal Node Forecasting: ADR-GNN\({}_{\text{T}}\)

The typical task in spatio-temporal datasets is to predict future quantities (e.g., driving speed) given several previous time steps (also called frames). Formally, one is given an input tensor \(\mathbf{X}_{temporal}\in\mathbb{R}^{n\times\tau_{in}c_{in}}\), where \(\tau_{in}\) is the number of input (observed) time frames, and the goal is to predict \(\tau_{out}\) time frames ahead, i.e., \(\mathbf{Y}_{\mathrm{temporal}}\in\mathbb{R}^{n\times\tau_{out}c_{out}}\). This is in contrast to'static' datasets such as Cora [59], where input node features \(\mathbf{X}\in\mathbb{R}^{n\times c_{in}}\) are given, and the goal is to fit to some ground-truth \(\mathbf{Y}\in\mathbb{R}^{n\times c_{out}}\). In this context, a'static' dataset can be thought of as setting \(\tau_{in}=\tau_{out}=1\) for the spatio-temporal settings. We show the overall flow of our 'temporal' architecture ADR-GNN\({}_{\text{T}}\) in Algorithm 45.

Footnote 5: In Algorithm 4, \(\oplus\) denotes channel-wise concatenation.

In our spatio-temporal ADR-GNN\({}_{\text{T}}\), we update the hidden state feature matrix \(\mathbf{U}^{(l)}_{\text{state}}\) based on the hidden historical feature matrix \(\mathbf{U}^{(l)}_{\text{hist}}\), as shown in Lines 6-9 in Algorithm 4.

Similarly to Attention models (see [C] in the additional appendix references), we incorporate time embedding based on the concatenation of sine and cosine function evaluations with varying frequencies multiplied by the time of the input frames, as input to our ADR-GNN\({}_{\text{T}}\), denoted by \(\mathbf{T}_{\mathrm{emb}}\in\mathbb{R}^{n\times\tau_{in}c_{t}}\), where we choose the number of frequencies to be 10, and by the concatenation of both sine and cosine lead to \(c_{t}=20\). We note that the time embedding is computed in a pre-processing fashion. To initialize the hidden feature matrices \(\mathbf{U}^{(0)}_{\text{state}},\ \mathbf{U}^{(0)}_{\text{hist}}\), we embed the input data \(\mathbf{X}_{\mathrm{temporal}}\), concatenated with \(\mathbf{T}_{\mathrm{emb}}\), using two fully connected layers, as described in Lines 3-4 in Algorithm 4. 6

Footnote 6: In Python notations, \(\mathbf{X}_{\mathrm{temporal}}[:,-c_{in}]\) extracts the last \(c_{in}\) entries of the second dimension of \(\mathbf{X}_{\mathrm{temporal}}\), which returns the features of the last time frame.

For the Chickenpox Hungary, PedalMe London, and Wikipedia Math datasets, we minimize the mean squared error (MSE) between the ground truth future node quantities and the predicted quanti
ties by \(\text{ADR-GNN}_{\text{T}}\), similar to the training procedure of the rest of the considered methods in Table 3. Specifically, following [68], the goal is to predict the node quantities of the next time frame given 4 previous time frames. On the METR-LA and PEMS-BAY datasets we minimize the mean absolute error (MAE), similar to [49], where we also follow the standard 12 previous time frames as inputs, and consider 3,6, and 12 future time frames node quantity prediction as output.

```
0: Node features \(\mathbf{X}_{\text{temporal}}\in\mathbb{R}^{n\times\tau_{in}c_{in}}\), time embedding \(\mathbf{T}_{\text{emb}}\in\mathbb{R}^{n\times\tau_{in}c_{t}}\) Output: Predicted future node quantities \(\tilde{\mathbf{Y}}\in\mathbb{R}^{n\times\tau_{out}c_{out}}\)
1:procedure\(\text{ADR-GNN}_{\text{T}}\)
2:\(\mathbf{X}_{\text{temporal}}\leftarrow\text{Dropout}(\mathbf{X}_{\text{ temporal}},p)\)
3:\(\mathbf{T}_{\text{emb}}\leftarrow\text{g}^{\text{time-embed}}(\mathbf{T}_{ \text{emb}})\)
4:\(\mathbf{U}_{\text{state}}^{(0)}=g_{\text{int}}^{\text{hist}}(\mathbf{X}_{ \text{temporal}}[:,-c_{in}]\oplus\mathbf{T}_{\text{emb}})\)
5:\(\mathbf{U}_{\text{hist}}^{(0)}=g_{\text{int}}^{\text{hist}}(\mathbf{X}_{ \text{temporal}}\oplus\mathbf{T}_{\text{emb}})\)
6:for\(I=0\ldots L-1\)do
7:\(\mathbf{U}_{\text{state}}^{(l)}\leftarrow\text{Dropout}(\mathbf{U}_{\text{ state}}^{(l)},p)\)
8: Advection: \(\mathbf{U}_{\text{state}}^{(l+1/3)}=\mathbf{U}_{\text{state}}^{(l)}+h\mathbf{ D}\mathbf{I}\mathbf{V}(\mathbf{V}(\mathbf{U}_{\text{hist}}^{(l)};\boldsymbol{ \theta}_{a}^{(l)})\mathbf{U}^{(l)}),\boldsymbol{\theta}_{a})\mathbf{U}_{ \text{state}}^{(l)}\)
9: Diffusion: \(\mathbf{U}_{\text{state}}^{(l+2/3)}=\text{mat}\left((\mathbf{I}+h\mathbf{K}^{(l )}\otimes\hat{\mathbf{L}})^{-1}\text{vec}(\mathbf{U}_{\text{state}}^{(l+1/3)})\right)\)
10: Reaction: \(\mathbf{U}_{\text{state}}^{(l+1)}=\mathbf{U}_{\text{state}}^{(l+2/3)}+hf( \mathbf{U}_{\text{hist}}^{(l+2/3)},\mathbf{U}_{\text{hist}}^{(0)};\boldsymbol{ \theta}_{r})\)
11:\(\mathbf{U}_{\text{hist}}^{(l+1)}=g_{l}^{\text{hist}}(\mathbf{U}_{\text{hist}}^ {(l)}\oplus\mathbf{U}_{\text{state}}^{(l+1)}\oplus\mathbf{T}_{\text{emb}})\)
12:endfor
13:\(\mathbf{U}_{\text{state}}^{(L)}\leftarrow\text{Dropout}(\mathbf{U}_{\text{ state}}^{(L)},p)\)
14:\(\tilde{\mathbf{Y}}=g_{\text{out}}^{\text{state}}(\mathbf{U}_{\text{state}}^{(L)})\)
15: Return \(\mathbf{Y}\)
16:endprocedure
```

**Algorithm 4**\(\text{ADR-GNN}_{\text{T}}\) Architecture Flow

## Appendix D Computational Complexity and Time

**Complexity.** Our ADR-GNN architectures include four main operations: (i) input/output embedding, (ii) advection, (iii) diffusion, and (iv) reaction layers.

The complexity of (i) and (iv) is \(\mathcal{O}(|\mathcal{V}|c^{2})\), where \(c\) is the number of channels, because they are composed of pointwise MLPs. The complexity of (ii) is \(\mathcal{O}((|\mathcal{V}|+|\mathcal{E}|)c^{2})\) because it requires the computation of the edge weights \(\mathbf{V}\), as shown in Algorithm 2, followed by a multiplication by the node features, as shown in Equation (5). Similarly, (iii) also is of complexity \(\mathcal{O}((|\mathcal{V}|+|\mathcal{E}|)c)\), because it is required to multiply the scaled Laplacian with the node features. Note that as discussed in Section 3.4, and similar to [14, 70], we do not explicitly invert the matrix \(\mathbf{I}+h\mathbf{K}^{(l)}\otimes\hat{\mathbf{L}}\), but rather use the conjugate-gradients (CG) method to solve a system of equations. Below, we further discuss the use of CG to solve the system of equations.

**Implicit solution of the diffusion term.** The diffusion step in our ADR-GNN, as shown in Section 3.4, requires the solution of the linear system at each step. As previously discussed this is solved by using the CG method. Thus, the backward (differentiation) function in most common software packages, such as PyTorch, tracks the CG iterations. This tracking can be avoidable by using implicit differentiation, which is the backbone of implicit methods (see [35] for detailed derivation). In the context of deep learning, implicit differentiation was used for implicit neural networks [33]. The basic idea is to use implicit differentiation of the equation

\[\mathbf{U}^{(l+2/3)}=\mathrm{mat}(((\mathbf{I}-h\mathbf{K}^{(l)}\otimes\hat{ \mathbf{L}})\mathrm{vec}(\mathbf{U}^{(l+1/3)})) \tag{13}\]

with respect to \(\mathbf{K}^{(l)}\) and thus avoid the tracking of the CG iterations if many are needed.

**Runtimes.** In addition to the complexity analysis above, we provide the measured runtimes in Table 6. As discussed in Section 5, learning the advection weights requires an increased computational effort. However, it can significantly improve the considered task metric. For convenience, in Table
6, in addition to the runtimes we also report the obtained task metric. Importantly, we show that the improved metrics offered by ADR-GNNs with 64 channels and 4 layers are not simply obtained due to the increased costs, by showing that enlarging GCN and GAT from standard 2 layers and 64 channels, to 2 layers and 256 channels (wide), or 64 layers and 64 channels (deep) does not yield similar improvements. We measure the runtimes using an Nvidia-RTX3090 with 24GB of memory, which is the same GPU used to conduct our experiments.

## Appendix E Hyperparameters

All hyperparameters were determined by grid search, and the ranges and sampling mechanism distributions are provided in Table 7. Note, that as discussed after Equation (8), we may add a BatchNorm layer before applying the non-linear activation \(\sigma\) to the reaction term, we therefore treat the use of batchnorm as a hyperparameter in Table 7.

## Appendix F Datasets

We report the statistics of the datasets used in our experiments in Table 8 and 9 for the node classification, and spatio-temporal node forecasting datasets, respectively. All datasets are publicly available, and appropriate references to the data sources are provided in the main paper.

\begin{table}
\begin{tabular}{l c c c c c c c} \hline \hline Metric & GCN & GAT & \begin{tabular}{c} GCN \\ (wide) \\ \end{tabular} & \begin{tabular}{c} GAT \\ (wide) \\ \end{tabular} & \begin{tabular}{c} GCN \\ (deep) \\ \end{tabular} & \begin{tabular}{c} GAT \\ (deep) \\ \end{tabular} & \begin{tabular}{c} GCN \\ (deep) \\ \end{tabular} & \begin{tabular}{c} GAT \\ (deep) \\ \end{tabular} & 
\begin{tabular}{c} ADR-GNNs \\ \end{tabular} \\ \hline Training time & 7.71 & 14.59 & 14.32 & 36.63 & 95.11 & 184.51 & 35.41 \\ Inference time & 1.75 & 2.98 & 2.86 & 7.57 & 12.93 & 38.96 & 8.24 \\ Parameters & 104 & 105 & 565 & 567 & 358 & 360 & 210 \\ Accuracy & 85.77 & 83.13 & 85.18 & 83.37 & 38.62 & 33.40 & 89.43 \\ \hline \hline \end{tabular}
\end{table}
Table 6: Training and inference GPU runtimes (milliseconds), number of parameters (thousands), and node classification accuracy (%) on Cora.

\begin{table}
\begin{tabular}{l c c} \hline \hline Hyperparameter & Range & Uniform Distribution \\ \hline input/output embedding learning rate & [1e-4, 1e-1] & log uniform \\ advection learning rate & [1e-4, 1e-1] & log uniform \\ diffusion learning rate & [1e-4, 1e-1] & log uniform \\ reaction learning rate & [1e-4, 1e-1] & log uniform \\ input/output embedding weight decay & [0, 1e-2] & uniform \\ advection weight decay & [0, 1e-2] & uniform \\ diffusion weight decay & [0, 1e-2] & uniform \\ reaction weight decay & [0, 1e-2] & uniform \\ input/output dropout & [0, 0.9] & uniform \\ hidden layer dropout & [0, 0.9] & uniform \\ use BatchNorm & \{ \textbackslash{}es / no \} & discrete uniform \\ step size h & [1e-3, 1] & uniform \\ layers & \{ \textbackslash{}2,4,8,16,32,64 \} & discrete uniform \\ channels & \{ 8,16,32,64,128,256 \} & discrete uniform \\ \hline \hline \end{tabular}
\end{table}
Table 7: Hyperparameter ranges
## Appendix G Experimental Results

### Additional Comparisons and Standard Deviations on Node Classification

To allow a more comprehensive comparison, and because some of the considered methods did not report the standard deviation around the mean accuracy, we now provide the experimental results from 4.1 on Cora, Citeseer, Pubmed datasets in Table 10, and Cornell, Texas, Wisconsin, Squirrel, Film, Chameleon datasets in Table 11,with the standard deviation around the mean of the 10 splits from [65]. Note, that here we do not color the tables, because some of the second or third top performing models in the main paper did not report the accuracy standard deviation, and therefore coloring Tables 10-11 will change the order of the best performing models.
also choose a random node \(\mathbf{v}_{\mathrm{dst}}\in\mathcal{V}_{\mathrm{ER}}\setminus\mathbf{v}_{ \mathrm{ER}}^{\mathrm{src}}\). The goal is to transport all the mass from all the nodes to the \(\mathbf{v}_{\mathrm{dst}}\), such that \(\mathbf{v}_{\mathrm{dst}}\) will have a feature of 1, the the rest of the nodes in the graph will be zeroes. That is, all the features in the graph are concentrated in \(\mathbf{v}_{\mathrm{dst}}\). In our example in Figure 2, we use the protocol specified here, to generate a graph with 5 nodes, and we show the approximation obtained with the advection, diffusion, and reaction terms.

### The Dirichlet Energy of ADR-GNN

We follow [71] and define the Dirichlet energy of the graph node features as:

\[\mathrm{E}(\mathbf{U}^{(l)})=\frac{1}{|\mathcal{V}|}\sum_{i\in\mathcal{V}} \sum_{j\in\mathcal{N}_{i}}||\mathbf{U}_{i}^{(l)}-\mathbf{U}_{j}^{(l)}||_{2}^{2}. \tag{14}\]

\begin{table}
\begin{tabular}{c c c c} \hline \hline Method & Cora & Citeseer & Pubmed \\ Homophily & 0.81 & 0.80 & 0.74 \\ \hline GCN & 85.77 \(\pm\) 1.27 & 73.68 \(\pm\) 1.36 & 88.13 \(\pm\) 0.50 \\ GAT & 86.37 \(\pm\) 0.48 & 74.32 \(\pm\) 1.23 & 87.62 \(\pm\) 1.10 \\ GCNII\({}^{\dagger}\) & 88.49 \(\pm\) 1.25 & 77.13 \(\pm\) 1.48 & 90.30 \(\pm\) 0.43 \\ Geom-GCN\({}^{\dagger}\) & 85.27 \(\pm\) 1.57 & 77.99 \(\pm\) 1.15 & 90.05 \(\pm\) 0.47 \\ MixHop & 87.61 \(\pm\) 2.03 & 76.26 \(\pm\) 2.95 & 85.31 \(\pm\) 2.29 \\ WRGAT & 88.20 \(\pm\) 2.26 & 76.81 \(\pm\) 1.89 & 88.52 \(\pm\) 0.92 \\ NSD\({}^{\dagger}\) & 87.14 \(\pm\) 1.13 & 77.14 \(\pm\) 1.57 & 89.49 \(\pm\) 0.40 \\ GGCN & 87.95 \(\pm\) 1.05 & 77.14 \(\pm\) 1.45 & 89.15 \(\pm\) 0.37 \\ H2GCN & 87.87 \(\pm\) 1.20 & 77.11 \(\pm\) 1.57 & 89.49 \(\pm\) 0.38 \\ LINKX & 84.64 \(\pm\) 1.13 & 73.19 \(\pm\) 0.99 & 87.86 \(\pm\) 0.77 \\ ACIII-GCN++ & 88.25 \(\pm\) 0.96 & 77.12 \(\pm\) 1.58 & 89.71 \(\pm\) 0.48 \\ \hline ADR-GNNs & 89.43 \(\pm\) 1.15 & 78.36 \(\pm\) 144 & 90.55 \(\pm\) 0.53 \\ \hline \hline \end{tabular}
\end{table}
Table 10: Node classification accuracy (\(\%\)) on _homophilic_ datasets. \(\dagger\) denotes the maximal accuracy of several proposed variants.

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline Method & Squirrel & Film & Cham. & Corn. & Texas & Wisc. \\ Homophily & 0.22 & 0.22 & 0.23 & 0.30 & 0.11 & 0.21 \\ \hline GCN & 23.96 \(\pm\) 2.01 & 26.86 \(\pm\) 1.10 & 28.18 \(\pm\) 2.24 & 52.70 \(\pm\) 5.30 & 52.16 \(\pm\) 5.16 & 48.92 \(\pm\) 3.06 \\ GAT & 30.03 \(\pm\) 1.55 & 28.45 \(\pm\) 0.89 & 42.93 \(\pm\) 2.50 & 54.32 \(\pm\) 5.05 & 58.38 \(\pm\) 6.63 & 49.41 \(\pm\) 4.09 \\ GCNII & 38.47 \(\pm\) 1.58 & 32.87 \(\pm\) 1.30 & 60.61 \(\pm\) 3.04 & 74.86 \(\pm\) 3.79 & 69.46 \(\pm\) 3.83 & 74.12 \(\pm\) 3.40 \\ Geom-GCN\({}^{\dagger}\) & 38.32 \(\pm\) 0.92 & 31.63 \(\pm\) 1.15 & 60.90 \(\pm\) 2.81 & 60.81 \(\pm\) 3.67 & 67.57 \(\pm\) 2.72 & 64.12 \(\pm\) 3.66 \\ MixHop & 43.80 \(\pm\) 1.48 & 32.22 \(\pm\) 2.34 & 60.50 \(\pm\) 2.53 & 73.51 \(\pm\) 6.34 & 77.84 \(\pm\) 7.73 & 75.88 \(\pm\) 4.90 \\ GRAND & 40.05 \(\pm\) 1.50 & 35.62 \(\pm\) 1.01 & 54.67 \(\pm\) 2.54 & 82.16 \(\pm\) 7.09 & 75.68 \(\pm\) 7.25 & 79.41 \(\pm\) 3.64 \\ NSD\({}^{\dagger}\) & 56.34 \(\pm\) 1.32 & 37.79 \(\pm\) 1.15 & 68.68 \(\pm\) 1.58 & 86.49 \(\pm\) 4.71 & 85.95 \(\pm\) 5.51 & 89.41 \(\pm\) 4.74 \\ WRGAT & 48.85 \(\pm\) 0.78 & 36.53 \(\pm\) 0.77 & 65.24 \(\pm\) 0.87 & 81.62 \(\pm\) 3.90 & 83.62 \(\pm\) 5.50 & 86.98 \(\pm\) 3.78 \\ MagNet & \multicolumn{3}{c}{ –} & 84.30 \(\pm\) 7.00 & 83.30 \(\pm\) 6.10 & 85.70 \(\pm\) 3.20 \\ GGCN & 55.17 \(\pm\) 1.58 & 37.81 \(\pm\) 1.56 & 71.14 \(\pm\) 1.84 & 85.68 \(\pm\) 6.63 & 84.86 \(\pm\) 4.55 & 86.86 \(\pm\) 3.29 \\ H2GCN & 36.48 \(\pm\) 1.86 & 35.70 \(\pm\) 1.00 & 60.11 \(\pm\) 1.71 & 82.70 \(\pm\) 5.28 & 84.86 \(\pm\) 7.23 & 87.65 \(\pm\) 4.98 \\ GraphCON\({}^{\dagger}\) & \multicolumn{3}{c}{ –} & \multicolumn{3}{c}{ –} & 84.30 \(\pm\) 4.80 & 85.40 \(\pm\) 4.20 & 87.80 \(\pm\) 3.30 \\ FAGCN & 42.59 \(\pm\) 0.69 & 34.87 \(\pm\) 1.35 & 55.22 \(\pm\) 2.11 & 79.19 \(\pm\) 5.87 & 82.43 \(\pm\) 2.87 & 82.94 \(\pm\) 1.58 \\ GPRGNN & 31.61 \(\pm\) 1.24 & 34.63 \(\pm\) 1.22 & 46.58 \(\pm\) 1.71 & 80.27 \(\pm\) 8.11 & 78.38 \(\pm\) 4.36 & 82.94 \(\pm\) 4.21 \\ ACMP-GCN & \multicolumn{3}{c}{ –} & \multicolumn{3}{c}{ –} & 85.40 \(\pm\) 7.00 & 86.20 \(\pm\) 3.00 & 86.10 \(\pm\) 4.00 \\ LINKX & 61.81 \(\pm\) 1.80 & 36.10 \(\pm\) 1.55 & 68.42 \(\pm\) 1.38 & 77.84 \(\pm\) 5.81 & 74.60 \(\pm\) 8.37 & 75.49 \(\pm\) 5.72 \\ GRAFF\({}^{\dagger}\) & 59.01 \(\pm\) 1.31 & 37.11 \(\pm\) 1.08 & 71.38 \(\pm\) 1.47 & 84.05 \(\pm\) 6.10 & 88.38 \(\pm\) 4.53 & 88.83 \(\pm\) 3.29 \\ G\({}^{2\dagger}\)\({}^{\dagger}\) & 64.26 \(\pm\) 2.38 & 37.30 \(\pm\) 1.01 & 71.40 \(\pm\) 2.38 & 87.30 \(\pm\) 4.84