Finite Volume Informed Graph Neural Network for Myocardial Perfusion Simulation

Authors: raoul  sallé de chou
matthew  sinclair
sabrina  lynch
nan  xiao
laurent  najman
irene e vignon-clementel
hugues  talbot
Venue: MIDL_io_2024_Conference
Published: 2024-02-01 12:26:47

Abstract: Medical imaging and numerical simulation of partial differential equations (PDEs) representing biophysical processes, have been combined in the past few decades to provide noninvasive diagnostic and treatment prediction tools for various diseases. Most approaches involve solving computationally expensive PDEs, which can hinder their effective deployment in clinical settings. To overcome this limitation, deep learning has emerged as a promising method to accelerate numerical solvers. One challenge persists however in the generalization abilities of these models, given the wide variety of patient morphologies. This study addresses this challenge by introducing a physics-informed graph neural network designed to solve Darcy equations for the simulation of myocardial perfusion. Leveraging a finite volume discretization of the equations as a "physics-informed" loss, our model was successfully trained and tested on a 3D synthetic dataset, namely meshes representing simplified myocardium shapes. Subsequent evaluation on a genuine myocardium mesh, extracted from patient Computed Tomography images, demonstrated promising results and generalized capabilities. Such a fast solver, within a differentiable learning framework will enable to tackle inverse problems based on $\text{H}_2$O-PET perfusion imaging data.

Main paper: 
 introduction: 
 In (Li et al., 2023) , a GNN with a finite volume constrained loss combined with ground truth data, successfully predicted the 2D Navier-Stokes solutions over a diverse range of shapes, starting from the solution of one time step and predicting the next one. (Donon et al., 2020) developed a FEM-informed loss function that effectively addressed the Poisson equation on various 2D shapes, employing a GNN without reliance on ground truth data. These neural networks demonstrated in some cases superior generalization abilities compared to datadriven-only models, but were only tested in 2D cases, with a small number of unknowns and with simple shapes.Our study introduces a Finite-Volume Informed Graph Neural Network for solving the 3D Darcy equations to model myocardial perfusion. It is trained on a synthetic dataset with diverse non-convex shapes resembling patient left ventricular geometries. Our model incorporates a physics-informed loss, thus eliminating the need for ground truth data. The PDE is discretized with Finite Volumes. Once the training is complete, the model can be applied in inference mode to new geometries to predict a solution, without retraining. We show preliminary results of our model and its generalization ability on a patient-specific leftventricle myocardium dataset from coronary computed tomography angiography (CCTA) images. To our knowledge, this is the first physics-informed model tested on a real-world 3D dataset with a clinical purpose, and the first fully-based finite volume-informed neural network (FVINN), hence without the need for ground-truth data. 
myocardial perfusion simulation: 
 In prior investigations, (Papamanolis et al., 2021 ) formulated a patient-specific simulation pipeline for diagnosing coronary artery diseases based on CT-scans. The model is initiated with the segmentation of the aorta, epicardial coronary arteries, and left ventricle myocardium from CCTA image data. To account for small arteries that fill the myocardium volume and are too small to segment on CCTA images, synthetic coronary trees are generated from the epicardial segmented roots (Jaquet et al., 2018) . Reduced order blood flow equations are then solved in the defined vascular network, with both segmented and synthetic components, to simulate the blood flow supply to the left ventricle myocardium.Subsequently, to simulate perfusion in the myocardium and to account for the microcirculation (blood flow in the arterioles, capillaries, . . .), a porous-media model governed by Darcy's law is utilized (Chapelle et al., 2010) . The model comprises a source term driven by coronary flow and a homogeneous sink term (to simulate an ideal venous system). In order to obtain the source term, p source , maps are generated by associating a perfusion territory in the myocardium with each outlet of the synthetic coronary tree. These territories are estimated through a Voronoi tessellation. Finally, a p source value is assigned to each perfusion region. This value is the pressure at the vasculature outlet, determined by solving the flow equations in the network.We will focus in this study on predicting the FVM solution of the Darcy equations on different meshes and p source maps, by a FVINN. We describe in the following subsection the different datasets for training, validation, and test.To facilitate model training, we generated a synthetic dataset comprising half-ellipsoidal shells aimed at simulating simplified myocardium morphologies. The dataset was created using gmsh (Geuzaine and Remacle, 2009) . As illustrated in Figure 1 , the geometries were generated by defining 5 points connected by elliptical arcs. The resulting 2D shape was then Figure 1 : Left: schematic of the synthetic shape generation. Centre: example of one final synthetic 3D mesh. Right: example of one patient myocardium with its p source map. extruded around the z-axis to form a volume, within which a volume mesh was generated. For the training dataset, 2000 distinct geometries were generated, with an additional 400 for the validation dataset. Various lengths L1, L2, L3, L4, and L5 were sampled to diversify the geometries. Furthermore, different p source maps were assigned to each geometry: within each volume, 1000 seed points were sampled. A Voronoi tessellation was grown from these points to emulate the perfusion regions associated with network outlets. Each Voronoi region was assigned a random p source value sampled from a normal distribution, in line with the observed p source distribution from patient-specific simulations (Papamanolis et al., 2021).As a test dataset, the left ventricle myocardium of 400 patients were segmented from CCTA image data, as part of the human-in-the-loop process for the HeartFlow FFR CT Analysis (HeartFlow Inc., US). CCTA images were received from multiple sites around the world by HeartFlow, Inc. (Taylor et al., 2013) . 3D tetrahedral meshes of the myocardium volumes were generated with gmsh. Synthetic p source maps were also attributed to each myocardium mesh according to the previously explained method 2.2. Each mesh was composed of approximately 60k nodes in both datasets. 
governing equations: 
 In this work we focus on solving the 3-dimensional (3D) single compartment Darcy model for an incompressible flow:w + K∇p = 0 in Ω (1) ∇.w = β source (p source -p) -β sink (p -p sink ) in Ω (2) ∂w ∂n = 0 in ∂Ω (3)with K the permeability tensor, w the Darcy velocity, p the capillary bed pressure, p source and p sink the source and sink pressure terms respectively, β source and β sink parameters describing the conductance of flow entering and exiting the myocardium respectively. These two parameters are patient specific; their computation is described in (Papamanolis et al., 2021) . By applying the divergence operator to Equation (1), we obtain the following Poisson equation, which can be solved for p with no-flux boundary conditions (Equation ( 3)):K∆p = -β source (p source -p) + β sink (p -p sink ) in Ω (4)FVM is a discretization method for the approximation of PDE. FVM is based on a volume integral formulation of the PDE over a finite number of control volumes discretizing the domain where the equations are solved in. The methodology unfolds by integrating the equations over each control volume. Introducing the vector n, normal to the control volume surface, and applying the divergence theorem on the diffusion term allows the replacement of volume integrals with surface integrals:EQUATIONIn this study, we consider the vertex-centered FVM where the control volumes are delimited by the dual mesh of a 3D tetrahedron element mesh (see Figure 2 ). Considering the mean value approach, the flux values (left-hand side of the equation here) are approximated by their values on the centre of each surface while the source (right-hand side) are approximated by their value at the centre of the control volume (p in the following is the discrete pressure):K f ∈faces ∇p f .n f S f = (-β source (p source -p) + β sink (p -p sink )) Ω (6)The values on each face are linearly interpolated between the nodes connected to the edge crossing the face f . First order approximations are used for the pressure gradient. GNNs represent a specific class of neural networks designed for processing graph structures. Specifically, meshes utilized for solving discretized PDEs can be viewed as graphs denoted as G = (N, V, E), where N is the number of nodes, V is a matrix of node features, and E is the adjacency matrix indicating connections between nodes. GNNs typically consist of multiple layers that update the node feature matrix by exchanging information between nodes along graph edges.The solution to the Darcy equations introduced in Section 3.1 necessitates solving a linear system of equations. Iterative methods, such as the Jacobi method, are very efficient to solve these large systems of equations. The Jacobi method starts with an initial guess of the solution and iteratively updates it based on the previous state:EQUATIONWith p k i the Darcy pressure of node i at iteration k, N i the indices of nodes belonging to node i neighbourhood, b i and a ij coefficients computed from the different constants in 6. The solution is obtained once the difference between two successive states is small enough. As mentioned by (Liu et al., 2021) , this sequence of operations is similar to the way GNNs perform.The graph attention network was first introduced by (Velickovic et al., 2017) . Given a graph G = (N, V, E), for every node i ∈ N , attention coefficients are computed only between neighbours: α ij = softmax j (LeakyReLU(a t [W h i ||W h j ]). With a a single-layer feedforward neural network, W a weight matrix and h i the vector of node features of nodes i. The attention mechanism finally updates the feature vector h i , similarly to the Jacobi method 7, through a nonlinear function σ:h i = σ( j∈N i α ij W h j ) (8)Similarly to (Zhao and Zhang, 2021), we employ a Graph Attention U-Net (GATU-net, Figure 3 ), with a slight deviation in the pooling and unpooling operations. The GATU-net architecture comprises three levels. During each pooling operation, a random sub-sample of the preceding upper-level mesh nodes serves as input. The node features of these subsamples are computed by averaging the node features of the k-nearest neighbours from the upper mesh's level. Notably, the unpooling function is the one proposed in PointNet++ (Qi et al., 2017) , where features are computed through distance-weighted interpolation from the k-nearest neighbours of the previous lower-level stage's mesh during the decoding phase.A node encoder is employed to encode node features into a latent space. For each node, the encoder takes as input the node coordinates, along with values for p source , β source , and β sink . Then, a node decoder processes the output of the GATU-net and predicts a pressure value for every node. The loss function to train the model is derived from Equation ( 7); it is the absolute mean over all nodes of the residual:|p θ (i) -( 1 a ii (b i - j∈N i a ij p θ (j)))| (9)With p θ (i) the pressure output of the node i of the decoder. It should be noted that no reference pressure solution is used to train the model. In addition, the same GATU-net architecture is trained with a classical L1 loss function between the network output and a ground truth from FVM. The latter model will be denoted as the supervised model (SM) in the next sections. The FVINN can be trained on a dataset comprising varying geometries and p source maps. The goal of the proposed method is to find a set of convolutions which minimize Equation ( 9) for every node of the meshes. Once trained, the model can be applied to a new geometry and p source map to predict the solution in inference mode. 
results: 
 In this section, we present the results of FVINN and SM compared to the FVM solution.For both models, we report the actual, relative and normalized mean absolute error (MAE) and Mean max absolute error per shape (MMAE). The normalized error is the actual error divided by the standard deviation of the FVM solution. The Computational time to solve the linear system of equations given in Equation ( 6) with using the GPU-enabled pytorch function linalg.solve is approximately 4.28 seconds. The forward pass of our network to compute the solution takes 0.56 second.The model was successfully trained on the synthetic dataset. SM obtained a slightly better MAE than FVINN on both the training and validation datasets. Specifically, on the validation dataset, the SM achieved a MAE of 3.74 Pa compared to 4.04 Pa for FVINN (Table 1 ). The relative MAE on the validation dataset was approximately 0.20% for both models, which corresponds to about 0.15 standard deviations of the FVM solution, as demonstrated by the normalized MMAE.On the test dataset, the MAE was also smaller for the SM compared to FVINN, with normalized MAE values of 0.40 and 0.33, respectively. However, FVINN exhibited better MMAE on all three datasets. Notably, the MMAE on the test dataset was 99.23 Pa and 35.10 Pa for SM and FVINN, respectively (Table 2 ). The high MMAE, combined with a relatively low MAE, indicates localized high errors for the SM, as displayed in Figure 4 . No particular locations of these errors were found.  
discussion: 
 In this study, we presented a graph neural network for 3D myocardial perfusion prediction.The model incorporates a finite volume-informed loss function and requires no ground truth data for training. Successfully trained on a fully synthetic dataset, the neural network was then tested on patient-specific myocardium meshes. To our knowledge, this is the first application of such a model on a dataset of patient-specific data and large size meshes. Despite our model not being trained against ground truth from simulations, it achieves comparable accuracies on synthetic datasets, validating this approach. However, while SM obtained better MAE, the MMAE of our FVINN was better on the three datasets and almost three times smaller than SM on the patient-specific dataset. Our model demonstrates more robust generalization results on unseen geometries, which may differ significantly from the original distribution of shapes in the training dataset. Nevertheless, the synthetic datasets used in our study represent highly simplified myocardium shapes, which limits our ability to draw conclusions about the efficacy of FVINN or SM for predicting Darcy solutions on patient left-ventricle myocardium meshes. Future work should involve training these models on patient-specific or more complex synthetic data to further investigate their performance.In conclusion, our study shows promising results for developing deep learning-based models as PDE solvers for real-world applications. Beyond reducing computational complexity, models such as the FVINN offer the potential to solve inverse problems from H 2 O-PET perfusion imaging data. Moreover, the integration of physics-based loss functions with an unsupervised framework holds potential for predicting more robust results while reducing the reliance on data, thereby supporting their application in medical context.  
 
 Appendix:  


Decision: This work innovates in designing a physics-guided GNN for 3D myocardial perfusion prediction. Thanks to the physics-constrained loss function, the method requires no ground truth during the learning process. This allows the model to be trained on a wide variety of p_source maps and geometries and therefore encourages its generalization to new cases. Result shows robustness of the model on unseen geometries and thus tends to validate the generalization hypothesis. The authors have successfully answered the reviewers' questions, which clearly improves the quality of their paper. 

For all these reasons, I have decided to accept this article.

Reviews: review: paper_summary: paper_summary: This paper presents a finite-volume based GNN to learn PINNs for Darcy equations to simulate myocardial perfusion.  The  PINN leverages the discretized formulation of  the Darcey equation over the tetreahedron mesh of the myocardium, as well as a GNN representing this geometry. The model were trained on synthetic datasets comprising simple half-ellipsoidal shells, and tested on myocardial meshes segmented from CCTA mage from patients.
main_review: main_review: The authors clarified a critical question for me from the original submission, regarding how the PINN can be trained for different geometry and source at the same time. This strengthens the novelty of the work. I'm thus increasing my rating.
strength_weakness: strength_weakness: The paper is missing some critical details about the methodology that prevents a proper assessment of the work proposed. See detailed comments below. 

Another primary concern of the paper is the lack of proper PINN baselines to assess the contribution for the work proposed. 

Numbers on computational cost for FVM vs. FVINN are necessary.
questions: questions: Please address the above concern about the methodology of the work (loss in (8)) and the training-testing setup.

Please add proper PINN baselines to demonstrate the benefit of 1) FV-informed formulations and 2) GNN formulations. If it is not possible to do a naive PINN on this, elaborate on why.

Please add computational cost for a FVM simulation vs. FVINN training to obtain the same 3D results.
paper_summary: paper_summary: The authors propose to train a physics-informed graph neural network (GNN) to solve the equations governing myocardial perfusion. The specific promise is that the model learns the physics purely from simulated data, yet generalises to real data sufficiently well. They show that this is qualitatively the case (the perfusion patterns match the ground truth obtained by "real" simulation), but quantitatively off by a factor.
main_review: main_review: The paper is well written, clear, and shows a combination of methods that are in my opinion a door opener for a variety of relevant simulation tasks in the medical field. Be it biophysical modeling of soft tissue deformation -- notoriously hard to parameterise properly and time consuming to actually calculate, or blood flow, or as here, perfusion: pointing out the need for research as this paper does for me is worth to have it presented to the community, hoping that it sparks some interest and discussion. 

Weak accept since I'm too little of an expert in this particular combination of methods, and in the underlying medical field to appreciate the contribution fully.
strength_weakness: strength_weakness: Bits and pieces I found unclear:
* The source map for the synthetic data is apparently completely random both in shape and in source pressure -- correct? If so, it might be stated. I was puzzled here, since I thought that this is more systematic in nature. I haven't worked in the field of cardiology, though, and can't judge this.
* The source maps were also generated randomly for the real data, so that only the shape of the myocardium seems to be different. Wasn't it possible to obtain data with a real segmentation of the coronary arteries and a better, more realistic estimate of their support regions?
* Modeling-wise: the vascular tree is not a mesh; so treating the entire myocardium like a sponge with arbitrary perfusion in all directions is perhaps not accurate; possibly it is good enough since the objective is just to simulate the last few millimeters from the ends of the vessels into the muscle fibers? My picture of the process, though, suggests that if a region with less supply is downstream of one with more supply, it will get more source pressure as the perfusion in the over-supplied region can't take up the supply -- the GNN on the other hand has only short-distance (next neighbor) interaction. Again, my mental picture might be wrong. I would like to hear a few words more on how the simplistic modeling approach of random source maps can be motivated clinically, though.
* It wasn't clear to me why the very few node features needed to go into the node encoder.
questions: questions: As I advocate acceptance of the paper as an oral, I don't have strong needs for a rebuttal.
paper_summary: paper_summary: The authors propose a novel method to accelerate the inference of myocardial perfusion simulations. Specifically, the method relies on a graph attention UNet, which is trained in an unsupervised manner by optimizing a physics-informed loss function derived from a finite-volume formulation. The method was trained on 2,000 synthetic  perfusion simulations with idealistic ventricular shapes and evaluated on additional simulations over 400 synthetic and 400 patient-specific geometries. The method achieves on both datasets a mean relative absolute error of 0.5% or lower.
main_review: main_review: I appreciate that the authors have addressed all of my questions and concerns. The changes made to the manuscript have greatly improved the quality of the paper. The authors have now added all necessary information to replicate the results. I am therefore increasing the rating.
strength_weakness: strength_weakness: The proposed method has several benefits, but I found several questions and concerns:

- First of all, the manuscript lacks essential information to replicate the results. Among other things, the paper is missing information on how the parameters for the synthetic geometries were sampled (e.g., range of parameters L1-L5 and the sampling pattern) as well as a clear description of the node encoder and decoder architecture is missing. Similarly, no hyperparameters for training the network were presented.

- Second, it seems to me like an essential limitation/requirement of the proposed method was not discussed at all. In particular, since the method takes as input nodal coordinates, the method would require the geometries to be in a canonical shape or it would require extensive augmentation to cover all possible rotations & translations. 

- Similarly, Figure 1 suggests that relatively homogeneous meshes were generated for these experiments. It is therefore unclear what impact mesh resolution has onto the prediction results. Could the authors please also comment on how the receptive field size of the network relate to the overall size of the underlying mesh?

- One particular claim of the proposed method is to have improved generalization capabilities thanks to the physics-informed/finite-volume-informed graph neural network. There is, however, no comparator. While I understand that a fair comparison is not straightforward, I strongly believe that it would be very valuable if the method is at least compared against a regular physics-informed neural network and/or a supervised graph deep learning approach. 

- In addition, the presented numerical evaluations only provide a global insight into the prediction accuracy. Neither the table nor the visual examples (since they are only showing a specific front of the geometry) can provide conclusions about any regional biases. Considering that cardiac ventricular geometries are considered, one simple next step may comprise the regional error analysis, e.g., based on the standardized 17-segment subdivision of the left ventricle by the American Heart Association. Moreover, it may be worth showing few more examples (front and back) of other heart geometries

- Finally, one particular advantage of such a method, as mentioned by the authors, may be the inference time once the network is trained. It would strengthen the paper, if a preliminary analysis on run-time differences between the ground truth finite-volume simulations and the network predictions (including the training time) were presented. In addition, I would appreciate if any information on the utilized software packages were presented.
questions: questions: In addition to the questions and concerns raised above, I would appreciate if the authors could comment on the following topics:

- The manuscript states in the conclusion "Moreover, the integration of physics-based loss functions with a non-supervised framework holds potential for enhancing the explainability of the models while reducing the reliance on data, thereby supporting their application in medical context". From my point of view, neither the unsupervised learning nor the physics-based loss helps in any way with explainability since there will be no guarantee that the output of the network is physically plausible. It may, however, help with reducing the errors, improved generalization, and faster convergence. Could the authors please share their opinion on this topic?

- Equation 3 presents the boundary condition of the Darcy model. Were any special considerations in the loss function taken to properly handle/incoporate this constraint?

- In the description of the graph pooling, it was mentioned that a random subset of nodes is chosen. Could the authors please clarify the reasoning behind it and how the connectivity for the coarse graph was computed?

Comments: title: Follow-up questions
comment: Thanks for the response. This makes more sense to me now. I'll update my rating.
FVINN can be trained on varying geometries and p_source maps at once as the node coordinates and p_source values are the input of our network. The goal of the proposed method is to find a set of convolutional operations that can reduce the loss function (8) for every node, knowing the position of node neighbors and their p_source values. Once we have learned the set of convolutional operations on the training set, we can apply it to new unseen geometries and p_source maps to predict a solution. The training and inference strategies are similar to classical supervised methods. In a way, given the mesh as input, the model needs to learn the coefficients a/b (which are computed using only finite element geometrical features and p_source values) so that it can successfully minimize the loss function during training and predict a solution at inference.
Thanks for your clarifications. Please elaborate how the proposed FVINN can be trained on varying geometries and p_source maps at once? In particular, how is the model trained as parameters a and b in the loss function (8) changes among samples, and how do new geometry or parameters a/b get incorporated to get a new PINN solutions on new geometry and p_source? These are critical innovations of the proposed method and are not clear from paper descriptions.
We thank the reviewer for her/his insightful comments and positive feedback. All the comments are addressed below in the same order as they appear in the review:

- The L1-L5 ranges were found empirically to represent the variability in patient LV volume ranges. The ranges in millimeters were [10, 20], [55, 85], [10, 20], [10, 20], [10, 20] for L1, L2, L3, L4 and L5, respectively. Uniform distribution within each range was used. The node encoder/decoder were simple feed-forward neural networks composed of two hidden layers with 64 neurons. Adam optimizer with initial learning rate of 0.001 and a 0.99 exponential decay was implemented. Batch of 8 meshes were passed to the network during training for 400 epochs. This information was added to the appendix.
- The meshes were all centered and aligned along the z-axis using, for each geometry, the axis passing by the apex and the shape barycenter. Even though this data normalization method is not perfect, it reduces the need for data augmentation to cover possible rotations & translations.
- We agree that the study of the impact of the mesh resolution should be conducted in further work. Experiments showed that approximately 30 graph convolutional blocks were necessary for the network receptive field to encompass the entire mesh. However, we found that, instead of using a single Graph Attention Neural Network with 30 convolutional blocks, a Graph U-net-like multiresolution architecture with 16 convolutional blocks on each level was performing even better. The pooling operations aim at increasing the receptive field. 
- We do not claim in the paper that we improve the generalization capabilities. We only talk about promising generalization capabilities: thanks to the FV-informed loss, the model shows good generalization on unseen geometry and p_source maps. This is in contrast to classical PINNs. We interpret those results as an improvement on the generalization capabilities of PINN. We realized an error on how we apply our activation functions in our model. Fixing the problem improved FVINN accuracy. The new results were added in the revised version of the paper. In addition, we provide a comparison of FVINN with a supervised model with the same network architecture implemented. While the supervised model obtained better MAE, it predicted higher MMAE overall on the all datasets. This last point indicates better robustness capabilities of our proposed FVINN compared to supervised models. 
- We agree that further work should comprise a regional error analysis based on the  standardized 17-segment subdivision of the left ventricle by the American Heart Association. Actually during the results analysis, many visual examples were analyzed and no regional bias was found. Additional visual results will be added in the appendix in the final version, if accepted. 
- It is well established that DL methods are faster at inference than classical numerical methods. To illustrate, we solve the linear system of equations in (7) using the GPU-enabled pytorch function linalg.solve the computational time is approximately 4.28s. The forward pass of our network to compute the solution takes only 0.56s. The model took approximately 2 days to train on two GPUs. An in-house code was used as the FVM solver. We have added this information in the revised version of the paper.

- Indeed, the Darcy velocity is a 3D vector. Modifications was done accordingly.
- Indeed, a minus was missing in equation (4).
- Changes so that both definitions are the same were made.
- The +- indeed denotes the standard deviation. Modifications were done accordingly.
- Modifications were done accordingly.
- Modifications were done accordingly.

- Thanks to the FVM-derived loss function, we can constrain more what is learning the network and therefore increase the potential for the model to predict physically plausible solutions.
- The boundary conditions mentioned in the Darcy model are incorporated in the model as source term on the mesh boundary element. As it is a zero-flux boundary condition, the source term to add is equal to zero. Therefore, no specific consideration in the loss function is required.
- A random subset of nodes from the less coarse levels were chosen for memory and computational efficiency. Similarly to classical numerical methods, we could have generated coarser meshes for each level of the GATU-net. However, this would induce more computational cost to generate multiple meshes. This method would also require to store for every data, all the mesh nodes coordinates, therefore increasing the GPU memory usage and decreasing the training batch size which would affect the performance of the network. The coarse mesh connectivity is computed as follows: clusters of nodes are generated by assigning every node of the upper-level to the closest node of the sub-sample. Then, edges are built between two sub-sample nodes if their clusters were connected by an edge in the upper-level.
We thank the reviewer for her/his insightful comments and positive feedback. We found a mistake in our code after the submission. Fixing the problem significantly improved the results. The modification was added to the revised version of the paper. All the comments are addressed below in the same order as they appear in the review:

- Indeed, the p_source maps and geometries are generated by random processes. We have clarified this point in the revised version of the paper. 

- We agree that the model should be tested on more realistic p_source maps. Further work is however needed to collect realistic p_source maps, e.g.  from the coronary artery segmentations and computational fluid dynamics simulations. We believe this issue is beyond the scope of the present paper.

- Indeed, the single compartment Darcy model should be accurate enough, as it aims at simulating only the blood flow in the small vessels such as arterioles and capillaries. A much more detailed vascular tree can be simulated from the end points of the segmented vessels (Jaquet et al.). Nonetheless, this perfusion model delivers accurate results, if it is coupled with blood flow simulation in the larger coronary arteries, using Navier-Stokes equations or reduced models of blood flow, as described in section 2.1. The coupled model requires iteratively solving flow in the vascular tree and in the myocardium, with different p_source maps obtained at each  iteration, until convergence, so that blood flow is well distributed between far-away regions. We refer the reviewer to the reference Papamanolis et al. for more information about the coupled model. Even though classical GNN only 
 have short-distance interaction, the Graph Attention U-net-like architecture we use allows for more distant interactions between nodes.

- The main advantage of the random p_source map is that they are quite easy to generate, and requires no additional simulation or segmentation, compared to more realistic ones. Those synthetic p_source maps were however generated within physiological ranges of p_source from Papamanolis et al. Their spatial heterogeneity is higher than the one derived from the coupled model described in the previous answer. Therefore, we believe that, if our model can predict accurate solutions on these synthetic maps, it should be able to generalize well on realistic ones.

- We agree with the reviewer that, with regard to the small number of node features, the node encoder could appear unnecessary. However, experiments do show that results are better with it. One possible explanation is that it could encode the relative position of the nodes to some specific region of the mesh ( this is plausible as the different meshes have overall similar shapes characteristics, and that b_source is proportional to the mesh volume).  We believe that it can help a node take into account whether or not it is close to a boundary, or to the apex, as the finite element, and therefore, the FV control volumes can be different.
We thank the reviewer for their insightful comments and positive feedback. We would like to inform the reviewer that a mistake was found in our code after the submission. Fixing the problem significantly improves the results. The changes was added and discussed in the revised version of the paper. We reply to the comments below in the same order as they appear in the review:

- Our FVINN can be trained on varying geometries and p_source maps both at once. This is not possible with classical PINNs. The trained FVINN model is then able to make predictions on new geometry with new boundary conditions at inference time. Indeed, similarly to supervised learning where the loss is computed using multiple ground truths as examples, the loss in (8) comprises patient-specific elements such as p_source values. The model learns a set of convolutional operations which can minimize the FV-informed loss starting from multiple p_source maps and geometries. The main difference between our method and classical supervised methods is that we use a more constrained loss function from FVM. Using these additional constraints, our model can predict a solution without any ground truth supervision.

- Similarly to the previous answer, we believe that our FVINN can not be directly compared to a PINN baseline. However, we added in the revised paper a comparison between FVINN and a classical supervised model. FVINN achieves comparable results with the supervised model, while not requiring any ground truth obtained from classical FEM/FVM solvers or actual measurements. Moreover, the results on the test dataset seem to indicate that FVINN has more robust generalization capabilities compared to the supervised method. The GNN formulation is implemented in this paper in order to predict Darcy solutions on unstructured grids like a FVM and infer on different geometries. To the best of our knowledge, GNN is the only common architecture which can be applied to unstructured grids. We refer the reviewer to the introduction of the paper for the explanations and resources on the advantages of GNN for our application.

- Indeed, our article lacked a computational cost comparison between the two methods. The computational time to solve the linear system of equations given in (7) with the GPU-enabled pytorch function linalg.solve is approximately 4.28 seconds. The forward pass of our network to compute the solution only takes 0.56 second. The two experiments were conducted on the same GPU. This comparison was added to the revised version of the paper.