Graph Fourier Neural Kernels (G-FuNK): Learning Solutions of Nonlinear Diffusive Parametric PDEs on Multiple Domains

Thank you so much for your support of our work. Below are our responses to weaknesses and questions. Weaknesses: In Table 1, for the three methods (G-FuNK, FNO, GNN), the number of parameters can differ by order of magnitude, making it challenging to evaluate the improvement in performance. Thank you for this observation, we would like to note that this was something we tried to work around but there were a few issues. First, for a FNO with a single Fourier layer, with a width of 50, there are already over 500k parameters. Additionally, for the GNN model, there was a bottle neck for the model to fit on the GPU while training and in addition we observed that making the GNN model larger only marginally improved its performance. The statement "Our method predicts entire trajectories in under 1 second, significantly outperforming traditional numerical methods. For example, cardiac EP simulations typically take at least 15 minutes on 12 CPU cores for one given set of initial conditions." could be more precise, particularly regarding the numerical solution being compared. It seems the cardiac EP simulations from numerical methods serve as the high-fidelity solutions that are used to generate the training data and evaluate the error. However, G-FuNK learns on reduced modes ($k_{max}$), which may lead to limited accuracy. It would be more informative if the numerical solutions are computed on coarser mesh that achieve similar of accuracy as G-FuNK, or using the reduced modes. This value should be more precise, the issue with giving a precise number is that the simulation time changes from patient to patient (i.e. bigger heart is a larger mesh and can take longer). A more precise answer is an average of 13.2 minutes ranging from 8.12 to 21.97 minutes. This information will be added to the manuscript. Questions: Given that the paper focuses on PDEs on graphs, it seems that the "Multipole Graph Neural Operator for Parametric Partial Differential Equations" (MGKN), which claims to be mesh-invariant, would be a more suitable baseline for comparison. While the original MGKN framework does not explicitly tackle changing geometry, it seems it can still be applied in this setting, since only the graph is used as an input. There are also many similarities between these two approaches: the Fourier transform and inverse Fourier transform in G-FuNK play a role similar to kernel convolutions in MGKN, which are computed using the multipole algorithm. Additionally, both methods employ some form of truncation to make computation more tractable (e.g., limiting modes or long-range interactions). I would appreciate further discussion or comparison of these two approaches. Thank you for this question. The MGKN relies upon message passing which uses an aggregation framework of the edges that was observed to inhibit a GNN to learn the effects of anisotropic diffusion on a traveling wave. We discuss this within the text on lines 132-134 and furthermore the GNN we presented as a comparison also uses a message passing framework, further presenting the need to develop a graph-based network that can better handle anisotropic data. The examples in that paper were only on a 1D problem and a 2D linear problem. We were not able to figure out a feasible way to adapt their code for 3D time-dependent dynamics. Additionally, with their method, it is difficult and sensitive problem to obtain block low rank multiscale matrices with nonideal partitions of the surface.

Thank you for reviewing our manuscript. We have address your comments below, Weaknesses: 1. The limitations should be addressed. 2. The so-called graph fourier transform is actually spectral method, which need the eigenvectors pre-calculated first. This procedure could make the method not useful for large scale problem. Computational efficiency and scalability should be reported, including offline computational for the eigenvectors and online computation time. Comparison with FNO and GNN in terms of efficiency is also absent. Thank you for this point, we agree that as the problem goes to a larger scale, the calculation of eigenvectors becomes more difficult. We report on line 212 that the computational complexity is of the order $O(k_{max}^2 n j)$ for computing the $k_{max}$ lowest eigenvectors of an undirected graph. Additionally, we are able to use a subsampled domain and predict the solution of a higher resolution domain which shows the advantage of using the global eigenvectors. The computation of the eigenvectors for our case takes on average 8.58 seconds for the 200 lowest modes for each of the 25 cardiac atrial geometries (approximately 20k nodes in each graph). We report the online computation time for the predictions to be less than 1 second for the cardiac ep examples, in contrast to the numerical simulation time of anywhere from 8-22 mins per example. All three networks were similar in efficiency in that both training time and inference were similar and a discussion will be added to the paper to note this. 3. The novelty of graph fourier transform is limited. This is a widely studyed area in graph neural network community. While we agree with the ubiquity of GFT, we propose that the novelty in our method comes from representing the PDE learning problem using the spectrum of the weighted Graph’s Laplacian for a more generalizable FNO that can also incorporate anisotropic fiber fields into the edges naturally (so that their effects on the solution field do not have to be directly learned by the network). Our comparisons demonstrate that for diffusive PDEs, using this approach instead of message passing GNN based methods is advantageous for including the fiber fields without increasing the number of input parameters in our network. In this regard, G-FuNK is very light weight. Questions: 1. Why do the authors use the neural ODE model? Any gain from the specific model? Yes, we use the neural ODE framework for a few reasons. First, learning $\frac{du}{dt}$ via a neural ODE allows for easily applying Dirichlet boundary conditions since before integrating the points could simply be set to a value. Second, in cardiac electrophysiology (and many other systems) it is common to apply external stimulus at different timepoints so being able to arbitrarily apply a “non-learned” external stimulus is of a great advantage. Third, not being constrained to a constant time step is an advantage since at different time points within a simulation a smaller time step might be necessary for better stability. While using a neural ODE framework is not necessary in all cases, we believe it provides numerous advantages compared to the standard methods. We can also use different time stepping methods (implicit schemes for improved stability, higher order methods when necessary).

Thank you for your review of our manuscript. We appreciate your input. Below we have responded to the comments in the review. -We will add to the manuscript the requirement that the domains $\Omega_{\alpha}$ must be diffeomorphic (we certainly do not need conformal equivalence and do not have conformal equivalence in any of the examples). Thank you for pointing out this missing detail. The domains do indeed need to be diffeomorphic. We do not provide a theoretical justification, based on the current state of the literature on operator learning theory, most results are for vanilla frameworks on fixed domains. There are rarely any theoretical guarantees for approximation error on model-specific variants of the canonical operator learning problem at least that we are aware of. In the simplest case, our method reduces to a standard Fourier Neural Operator for which there are universal approximation guarantees by Kovachki et al., but we can of course handle more general domains than just a regular square grid. -In the case of the heart surface, the PDE is being solved by finite elements on the surface and the anisotropic Laplacian term is defined on the surface, which is precisely what we aim to approximate with our graph Laplacian on the surface. The main point of G-FuNK is that the canonical Fourier construction is replaced by the graph Laplacian which adapts to both the diffusivity fields and the surface on which they are defined. The specific architecture is an analog to the Fourier Neural Operators. -The PDE solver for the 3D Cardiac EP example, which is an open-source fully optimized finite element library for these problems, takes roughly 13 minutes to execute for a single set of initial conditions. For applications of these models, one is interested in doing parameter sweeps across initial conditions and stimulus to find some optimal treatment procedures for patients. It is unclear to us what the reviewer means by these being low dimensional problems. To our knowledge, the literature on neural operators for PDEs seems to be primarily focused on 2D and 3D examples, which is what we provide benchmarks for.

We would like to thank the reviewer for their comments and their time in reviewing our manuscript. Response to your questions are below. -For the heat equation example we used the software Mathematica which takes as an input the continuous domain, PDE and Initial/Boundary condition and solves the PDE to arbitrarily high precision. The rate of convergence is not applicable here since we do not pass in a discretization of the domain. (The numerical convergence examples were previously added to the Appendix) -There are several different errors which accumulate to give the final error (e.g., number of training samples, number of nodes in the graph, optimization error). The approximation of the eigenfunctions via the graph Laplacian is one of these errors and it is not likely that the analytical eigenfunction would improve the error more than the changes to the other factors mentioned above. In addition, our diffusion coefficient is not constant. Therefore, we don’t think it is feasible to compute the analytical eigenfunctions for each sample and since the purpose of this paper is to predict on anisotropic domains, we do not see the utility in these computations for our purposes.

We would like to thank the reviewer for their time and effort put into reviewing our manuscript. We addressed their comments. Major: W1) Some works provide theoretical approximation error for vanilla neural operator frameworks on a single domain, however, they do not provide generalization error bounds at all or explain the numerical performance to address. Even regarding theoretical results on approximation error, these generally are quite broad in terms of the class of models that permit and do not have specific rates. Model specific approximation error results are hard to come by, especially on problems defined on a family of domains. W2 and W6) Yes, incorporating a PDE loss, as in PINNs, is possible within this framework, but we chose not to include it to highlight the specific advantages of the G-FuNK layer. Adding a PDE loss would make it challenging to attribute the model's success to our method's novelties and significantly increase training time. For the cardiac EP example, implementing PINNs would be particularly complex due to the need for a local metric tensor or coordinate system on the 3D surface, as Euclidean space alone cannot represent the spatial derivatives. Additionally, parameterizing the multiply connected (N=5) topological domain of the cardiac geometry in $\mathbb{R}^3$ is non-trivial. Our approach is designed as a data-driven framework that can generalize across a family of PDEs without requiring explicit knowledge of the underlying equations. While our training data comes from simulations of known PDEs, the architecture leverages this indirectly, such as by incorporating diffusivity into the weighted graph. In the cardiac EP example, solution errors primarily arise from wavefront time delays, but the overall wavefront form remains physically consistent. PINNs are better suited for cases where physical principles cannot be learned from the data alone, but they may not be ideal for this application due to the challenges outlined above and their computational expense. W3) No, the solution does not need to be smooth. We present in both the reaction-diffusion and cardiac EP examples a solution where the “wavefront” is almost discontinuous and the PDE is nonlinear. Examples of this discontinuity are shown in Courtemanche et. al. It is expected that G-FuNK would perform similar to FNO/GEO-FNO in any examples such as the viscous Burgers equation previously shown in other works but with the advantage noted within this paper. We note that G-FuNK can be constructed to be identical to FNO but can also handle other cases. W4) This is a greatly discussed problem with neural ODEs and is an ongoing field of research. It should be noted that while the proposed work uses a neural ODE, G-FuNK can operate to predict the next time step like any other previous works have done. In particular, G-FuNK can be combined with an implicit solver, of course with significantly increased computational cost. W5) We appreciate the reviewer’s suggestion to compare data generation, training, and prediction times. The primary motivation behind G-FuNK is to accelerate clinical decision-making for new patient-specific geometries of the left atrium, where multiple pacing sites and frequencies must be tested to determine optimal ablation strategies. In a standard finite element solver workflow, each simulation for a patient-specific geometry takes approximately 13.5 minutes per pacing site. Testing multiple pacing sites for numerous patients quickly becomes computationally expensive. In contrast, G-FuNK requires a one-time offline training cost (1–2 days) on 24 patient geometries with varying pacing sites. Once trained, it can make predictions for new geometries and pacing configurations in seconds, enabling instantaneous parameter sweeps to determine optimal treatment strategies. This one-time effort amortizes over future test cases, significantly reducing the overall computational burden as the model generalizes without requiring retraining. We will clarify this in the revised manuscript and include a comparison of the computational costs in a clinical workflow, highlighting how G-FuNK achieves cost-effectiveness for large-scale simulation predictions, particularly for cardiac EP. W7) In example 2 and 3, we extrapolate to different domains and in example 1 and 3 we extrapolate to different parameters of the PDE. In fact, example 2 is almost exactly what you suggest with the different domains, but with varying sized 2D rectangles. The test domains for the reaction diffusion and cardiac EP examples are unseen in the training. The cardiac EP example shows a case of extrapolation and domain-adaptation as the test domain is on a never-before-seen geometry with completely different boundaries and domain. We also believe this example shows that it is parameter adaptable since the test geometry has its own unique diffusive field that is different from all the other domains used in training. Will fix minor issues.