Deep learning of material transport in complex neurite networks

Neurons exhibit complex geometry in their branched networks of neurites which is essential to the function of individual neuron but also brings challenges to transport a wide variety of essential materials throughout their neurite networks for their survival and function. While numerical methods like isogeometric analysis (IGA) have been used for modeling the material transport process via solving partial differential equations (PDEs), they require long computation time and huge computation resources to ensure accurate geometry representation and solution, thus limit their biomedical application. Here we present a graph neural network (GNN)-based deep learning model to learn the IGA-based material transport simulation and provide fast material concentration prediction within neurite networks of any topology. Given input boundary conditions and geometry configurations, the well-trained model can predict the dynamical concentration change during the transport process with an average error less than 10% and 120~330 times faster compared to IGA simulations. The effectiveness of the proposed model is demonstrated within several complex neurite networks.


Introduction
The geometry of neurites is known to exhibit complex morphology which is essential for neuronal function and biochemical signal transmission. However, the highly branched networks of neurites also make it challenging to mediate the intracellular material transport because the material synthesis and degradation in neurons are carried out mainly in the cell body 1, 2 which leads to a long-distance transport for the material. The disruption of this long-distance transport can induce neurological and neurodegenerative diseases like Huntington's, Parkinson's and Alzheimer's disease [3][4][5][6][7] . Therefore, it has attracted considerable amount of attention in recent years to study the transport mechanisms and build mathematical transport models in neurite networks. Recent studies show that the molecular motors play fundamental roles in the intracellular transport to carry the material and move directionally along the cytoskeletal structure like microtubules or actin filaments [8][9][10] .
Motivated by these findings, different mathematical models based on partial differential equations (PDEs) have been proposed to help understand the transport mechanisms and the pathology of neuron diseases. For instance, Smith and Simmons developed a generic model of the molecular motor-assisted transport of cell organelles and vesicles along filaments 11 . Based on this motor-assisted transport model, Friedman and Craciun presented an axonal transport model by considering the ability of material to bind more than one motor protein 12 . Craciun et al. introduced the pausing of the material during the transport process in the model to study the slow axonal transport 13 . In addition, several PDE models were developed to study the transport impairment by accounting for the traffic jam 14 and microtubule swirls 15 in abnormal neurons. Though the aforementioned models provide reasonable mechanistic explanation of the transport mechanism or the formation of the transport impairment in neurites, most of these models were solved only in one-dimensional (1D) domain without considering the impact of complex neurite geometry. Recent advances in numerical methods and computing resources allow us to simulate detailed cellular process in complex geometry using 1D or three-dimensional (3D) PDE models. For instance, computational softwares based on finite element method (FEM) have been used to solve PDE models in neuron 16 and cell biological systems 17 . However, these tools have accuracy and computational cost issues when tackling highly branched geometry like neurite networks. Based on the conventional FEM, isogeometric analysis (IGA) 18 was proposed to directly integrate geometric modeling with numerical simulation and achieve better accuracy and robustness compared to FEM. With the advances in IGA, one can simulate the transport process by solving the PDEs in an accurately reconstructed neuron geometry. In our previous study, we developed an IGA-based simulation platform to solve a 3D motor-assisted transport model in several complex neurite networks and obtain the spatiotemporal material distribution 19 . However, the high computational cost to perform 3D simulations has limited its application in the biomedical field when fast feedback from the computer simulation is needed. Figure 1. An overview of the GNN model to learn and predict neuron material transport process. The input neurite network is first decomposed into pipes and bifurcations to create the graph representation of the neurite network. Next, the input features x i of each pipe or bifurcation are processed by the corresponding GNN simulator (G p S or G b S ) to generate the intermediate concentration result x mid . Then, the GNN assembly model (G A ) takes x mid as input and compute the interaction between simulators to predict concentration result x o .
To address the limitation in the current simulation platform, we propose to build a surrogate model by combining deep learning (DL) with IGA simulation. DL has been proven successful in computer vision and image recognition [20][21][22] by handling large amount of labeled datasets and providing fast end-to-end prediction. The practical success of DL in artificial intelligence also inspires its application in solving high-dimensional PDEs 23 and learning the physics behind PDE models 24 . In particular, deep neural networks (DNNs) are becoming popular in surrogate modeling because it can be far more efficient when predicting complex phenomena 25 . For instance, Farimani et al. applied the conditional generative adversarial network (cGAN) in a data-driven paradigm for rapid inference, modeling and simulation of transport phenomena 26 . Wiewel et al. proposed a DL-based fluid simulator to predict the changes of pressure fields over time 27 . Li et al. developed an encoder-decoder based convolutional neural network (CNN) to directly predict the concentraion distribution of a reaction-diffusion system, bypassing the expensive FEM calcuation process 28 . While these works manage to learn the underlying physical models for prediction, they are limited to handle the problem in relatively simple geometry with Euclidean data (e.g. structured grid) available for training. Recently, many studies on graph neural networks (GNNs) have emerged to extend DL techniques for data defined in graph structure. Inspired by the success of CNNs 20 , a large number of methods were developed to re-define the convolutional operation for graph data and achieve great performance in computer vision tasks like graph node classification and image graph classification [29][30][31][32] . GNNs were also applied to predict the drag force associated with laminar flow around airfoils 33 or the property of crystalline materials 34 , understand the interaction behind physics scenes 35 , solve PDEs by modeling spatial transformation functions 36 , or learn particle-based simulation in complex physical systems 37 .
In this study, we develop a GNN-based model to learn the material transport mechanism from simulation data and provide fast prediction of the transport process within complex geometry of neurite networks. The use of GNN is motivated by the extensive topologies of neurite networks and the IGA simulation data stored in mesh structure. To ensure the model is applicable to any neurite geometry, we build the graph representation of the neuron by decomposing the neuron geometry into two basic structures: pipe and bifurcation. Different GNN simulators are designed for these two basic structures to predict the spatiotemporal concentration distribution given input simulation parameters and boundary conditions. Specifically, we add the residual terms from PDEs to instruct the model to learn the physics behind the simulation data. To reconstruct the neurite network, a GNN-based assembly model is used to combine all the pipes and bifurcations following the graph representation. The loss function of the assembly model is designed to impose consistent concentration results on the interface between pipe and bifurcation. The well trained GNN model can provide fast and accurate prediction of the material concentration distribution, leading to an efficient DL framework for studying the transport process in complex 3D models of neurite network. The framework was applied to several complex neurite networks achieving an average error less than 10% and 120 ∼ 330 times faster compared to IGA simulations.

GNN model overview
The aim of our GNN model is to learn from the material transport simulation data and predict the transport process in any neurite network. However, the geometry diversity of the neurite networks makes it impossible to train the DL model directly on the entire neurite network. To address this issue, we introduce the graph representation of the neurite network and build the GNN model based on the graph network (GN) framework 35 . A neurite network can be decomposed into two basic structures: pipe and bifurcation (yellow square and red circle in Fig. 1, respectively). Each structure can be treated as one node in the graph and the nodes can be connected following the skeleton of the neurite network to constitute the graph. Based on the graph representation of the neurite network, two separate GNN simulators are built for the pipe (G p S ) and the bifurcation (G b S ), respectively. Given input features x i including node locations, simulation parameters and initial nodal concentration, the simulators can output the intermediate nodal concentration result x mid in the pipe and the bifurcation, respectively. To obtain a consistent global concentration result, a GNN assembly model (G A ) is used to learn the interaction between different structures so that given intermediate value We implement our GNN model using PyTorch 38 and the "PyTorch Geometric" library 39 . The detailed results will be explained in the following sections.

GNN simulators for pipe and bifurcation
Since the pipe and the bifurcation have different geometry topologies, we train two separate simulators to handle different graph structures extracted from the simulation results. Both simulators share the same recurrent "GN block + MLP Decoder" architecture but are trained with pipe and bifurcation datasets, respectively. As shown in Fig. 2, the input feature vectors depicting geometry, parameter settings, boundary conditions and concentration values c t k at t k are first processed by a series of GN blocks to learn the hidden layer representations that encode both local graph structure and features of nodes. In our GNN simulator, the GN block includes two modules φ e , φ v to update the edge and node features, respectively. φ e computes the concentration gradient along the edge and the length of the edge. φ v is a multilayer perceptron (MLP) consisting of two hidden layers (with ReLU activations), each layer with the size of 32. The forward computation of our GN block is characterized by Algorithm 1. After the GN blocks output the latent graph that encodes all the node feature, a MLP is then used to decode the hidden layer representation and output predicted concentration values c t k+1 at the next time step t k+1 . The MLP has three hidden layers (with ReLU activations), followed by a non-activated output layer, each layer with the size of 32.

Algorithm 1: GN Block
Input: for each edge e i j ∈ E do Compute edge attributes: Aggregate edge attributes: x e,v i = ∑ j x e i j ; Compute output: The input graph of the simulator is created by extracting certain nodes from each cross section of a hexahedral mesh following the templates ( Supplementary Fig. S1). For each pipe structure, we use the circular plane template ( Supplementary  Fig. S1C) to extract 17 nodes from each cross section. For each bifurcation structure, we use another template ( Supplementary  Fig. S1D) to extract 23 nodes from the cross section at the branch point of the bifurcation. We also use the same circular plane template to extract 3 circular cross sections on each branch of the bifurcation. For each node, we collect the geometry information and simulation parameters in an input feature vector. Here, the geometry information of each node is encoded by its coordinates and the radius of the cross section on which the node is located. The nodal concentration value is set to be the target prediction.
To encourage the model to learn the physics underlying the simulation data, we add the residuals of the governing equation where c 0 , c + and c − are the spatial concentrations; D is the diffusion coefficient of materials; u u u + and u u u − are velocities of materials; k ± and k ± are filament attachment and detachment rates (See Eq. 3 in Methods for the detailed physical meaning of the aforementioned variables); N denotes the number of nodes on the graph; and superscripts P and G denote the prediction and the ground truth value on the graph, respectively. To create the training dataset for two simulators, we first use our IGA solver to run the material transport simulation (see Methods for details) in two complex zebrafish neurons from the NeuroMorpho.Org database 40 (NMO_66731 and NMO_66748 in Fig. 3B&D). Regarding the simulation setting, we use 200 different boundary condition values and set constant parameters in Eq. 3 and Eq. 5 as D = 1.0 µm 2 /s, k ± = 1.0 s −1 , k ± = 0.5 s −1 and u i = 0.1 µm/s. The time step for each simulation is set to be 0.1 s. We simulate until the transport process is steady and then extract the spatiotemporal simulation results of 100 different geometries for each simulator from these two neurons. The nodal feature vectors and nodal concentration values are stored for each geometry to build the training dataset that contains 20,000 samples for each simulator. After establishing the architecture of two simulators, we train each simulator by randomly selecting 75% samples as the training data. The Adam optimizer 41 is used to optimize the model with the step learning rate decay ranging from 10 −3 to 10 −6 . Our simulators are trained for 200 4/18 epochs and we evaluate the performance of the model using the rest 25% samples as the test dataset. At the end of training, the test loss for the pipe and bifurcation simulators converges to 0.265 and 0.111, respectively ( Supplementary Fig. S2). To demonstrate the performance of our GNN simulators, we select three prediction results of a pipe and a bifurcation from the test dataset and compare with their corresponding IGA simulation results in Fig. 2. For the pipe simulator, we find the model can accurately capture the boundary conditions and predict the concentration distribution at each time step, which indicates that the GNN simulator manages to learn the time-dependent behaviour of the transport equations. By comparing nodal error at t = 2.6 s and t = 3.4 s in Fig. 2B, we find the error increases along with the front of material propagation through the pipe, which suggests that the pipe simulator is sensitive to the sudden distribution change during the material propagation and needs further improvement. For the bifurcation simulator, the predicted result has higher accuracy around the branch region which suggests that the bifurcation simulator learns to transport the material in correct directions at the branch point. However, the boundary condition is not preserved as well as the pipe simulator performs. The possible reason is that the boundary nodes have less neighbouring nodes for edge feature aggregation compared to the interior nodes and the lack of neighbouring information 5/18 leads to higher error. Therefore, our bifurcation simulator can be further improved to balance the neighbourhood between the boundary and interior nodes to better preserve the boundary condition.
We perform 4-fold cross validation and the mean relative error (MRE) results for both simulators are shown in the second column of Table 1. The well-trained GNN simulators can provide concentration prediction with an error of less than 8% in average. With the standard deviation less than 1.3%, the model shows good performance when being generalized to unknown data. To study the impact of the PDE residuals in Eq. 1 on the prediction performance, we compare our simulators with the model trained using the standard MSE loss function. We perform cross validation with the same dataset for each comparison and the results are shown in Table 1. The comparison shows that the model achieved better accuracy when trained with the "MSE + PDE residuals" loss function, which indicates the physics information contained in PDE residuals is learned by the model and improves its performance.

GNN assembly model
The objective of the GNN assembly model is to assemble the local prediction from simulators and output an improved continuous concentration distribution on the entire geometry. An overview of the GNN assembly model is shown in Fig. 3. The assembly model needs to learn the interaction between the simulators. Here the assembly model includes three components to cover all the assembly scenario in a decomposed neurite network: (a) pipe and pipe G p−p A ; (b) pipe and bifurcation G p−b A ; and (c) bifurcation and bifurcation G b−b A . During the assembly, the model loops each simulator node on the decomposed neurite network and utilizes the "message-passing" scheme to gather the predicted results from its neighbouring simulator nodes (green arrows in Fig. 3A). In particular, the nodal predicted results on the interface between two simulator nodes are collected. Then, all the collected values from the neighbouring simulator nodes are concatenated with the values from the current simulator node and processed by a MLP to improve the prediction. While have different numbers of input and output values, they share the same MLP architecture with three hidden layers (with ReLU activations), followed by a non-activated output layer, each layer with the size of 32.
To ensure the concentration result is consistent on the interface between two simulators, we add a penalty term to the MSE loss function as where N denotes the number of nodes on two assembled simulators, M denotes the number of nodes on the interface, superscript P and G denote the prediction and the ground truth value on the graph, respectively. Superscripts s 1 and s 2 denote the prediction value from the first and second simulators, respectively. α is the penalty strength which is set to be 10 in this study. We use the same IGA simulation results in two complex geometries of zebrafish neurons (NMO_66731 and NMO_66748 in Fig. 3B&D) to create training dataset for the assembly model. Based on the graph representation of these two trees, we follow three basic assembly structures and extract 30 different geometries for each type of assembly. The simulation results of these geometries are output and create 6,000 samples for each assembly structure. The model is trained using 75% samples as the training dataset and we evaluate the performance of the model using the rest 25% samples as the test dataset. At the end of training, the test loss for the GNN assembly model converges to 0.1492 (Supplementary Fig. S3). To test the performance of the assembly model, we use these two complex geometries and compare the concentration prediction at steady state ( Fig.  3B-E). The prediction MREs are 6.7% for NMO_66731 and 7.3% for NMO_66748, respectively. We find that the assembly model manages to generate a continuous global concentration distribution for each geometry. The comparison of predicted concentration between nodal error for each geometry also shows that the high concentration regions are predicted accurately. This indicates that our GNN model can capture the locations that are easier to develop transport disruption due to material accumulation.

Results for complex neurite networks
After the simulators and the assembly model in the GNN framework are well trained, we use our GNN model to predict the concentration distribution in several complex neurite networks and compare with the simulation results from our IGA solver. All complex neurite networks selected from the NeuroMorpho.Org database are shown in Figs. 4 and 5. Since the GNN model is trained with the simulation data in zebrafish neurons, we pick another two zebrafish neurons (Fig. 4A&B) to test the model performance in the neurons from the same species. We also pick another six mouse neurons (Fig. 4C, 4G, 4H and 5) from a different species to validate the model. We choose neurons with the number of bifurcations ranging from 9 to 356. For each neuron, we first run IGA simulation to get the ground truth concentration results and then compare with the predicted concentration obtained from our GNN model. Fig. 4 shows the prediction results for the neurons with longer branches and Fig. 5, shows neurons with more bifurcations and shorter branches. We observe that our GNN model predicts the distribution well with relative larger error around the inlet regions with high concentration. In addition, the model performance gets worse 7/18 in longer branches or regions with high density of bifurcations. The possible reason is that the increasing complexity of the geometry leads to a larger graph representation which affects the prediction accuracy due to the complicated assembling process on more simulators. The detailed geometry information, computation statistics and error comparison are summarized in Table 2. By comparing the prediction MRE of zebrafish neurons (top four rows in Table 2), we find the MRE only increases by around 0.9% for all the tested zebrafish neurons, which validates the applicability of the model among neurons from the same species.
In addition, we find that the average prediction MREs of zebrafish and mouse neurons are comparable with 7.25% and 8.5%, respectively, indicating that the trained model can also work for neurons from other species. To evaluate the computation performance, we define a speedup ratio as the ratio between the IGA computation time and the GNN prediction time. We run the code and measure the computation time on the XSEDE (Extreme Science and Engineering Discovery Environment) supercomputer Bridges 42 in the Pittsburgh Supercomputer Center. For each geometry, we perform the IGA simulation through parallel computing on CPU and perform GNN prediction on GPU. The speedup ratio for each geometry is shown in Table 2 and we find that our GNN model can achieve up to 330 times faster compared to IGA simulation. We also observe that the speedup ratio decreases when the neuron geometry becomes complicated with more bifurcations. The ratio converges to around 120 when the bifurcation number reaches 356, which is still a significant improvement by reducing computation time from hours to minutes.

Discussion
In this paper, we present a GNN-based DL model for predicting the material transport process in complex geometry of neurons. Our strategy utilizes a graph representation of neurons to overcome the difficulty in handling different neuron geometries with the DL model. The underlying assumption of our approach is that the material concentration distribution can be predicted locally in simple geometries and then be assembled following the graph representation to restore the concentration distribution 8/18  (Fig. 5E) (1,350,864, 1,179 in any complex geometry. Two GNN models are developed to validate our assumption. The first GNN-based model serves as the simulator to predict the dynamic concentration results locally on the mesh of two basic structures: pipe and bifurcation. We adopt GNN to directly handle the IGA simulation data stored in the unstructured mesh. Given input boundary conditions and geometry information of a pipe or bifurcation, the well trained simulators are able to provide high accuracy results (MRE < 7%). The second GNN-based assembly model then collects all the local predictions and follows the graph representation of the neurite network to update the global prediction. As shown in Figs. 4-5, the model is evaluated on complex neurons from mouse and zebrafish and shows its applicability with consistent performance in different neuron geometries. In particular, the model is capable of providing the spatiotemporal concentration prediction with MRE < 10% and over 100 times faster than the IGA simulation. Furthermore, the accurate distribution prediction can help us determine high concentration regions, which is essential to infer possible locations to develop transport disruption.
Our study shows that the complex and diverse geometry of neurons has major impact on the material concentration distribution and further affects the prediction accuracy of our GNN model. As shown in Figs. 4-5, the radius of most complex neurons is larger around the inlet region and decreases downstream, which contributes to material accumulation at the inlet region. The increase of bifurcations can also aggravate the accumulation when there is sharp decrease of radius in the downstream branches (Fig. 5E). These geometry features contribute to the complicated concentration distribution and thus bring challenges to improve the performance of our GNN model. One challenge we have for the most complex neuron NMO_00865 is the MRE gets worse to 12.5% when it was initially tested using the model trained with zebrafish neuron dataset. We plot the nodal error and find the error is higher at regions with high curvature or sharp radius change (Supplementary Fig.  S5). The possible reason is the GNN model lacks the knowledge of these scenarios since they are not common in the training dataset obtained from zebrafish neurons. To address this issue, we extract 20 geometries for each scenario from NMO_00865 and include them into the training dataset. We adopt the transfer learning method 43 to reuse the pre-trained GNN model as the starting point and obtain an improved model by training with the new dataset. The new model achieves the MRE of 9.2% on NMO_00865 which is comparable to the other testing neurons. This indicates that our GNN model can be further optimized with transfer learning on a better training dataset and an optimal dataset with a variety of geometries considered is quite essential to improve the applicability of the model.
The integration of the governing equations with the GNN model also plays a critical role in guiding the GNN model to learn the underlying physics behind the simulation data. By including the PDE residuals as the physics loss in the loss function of the GNN simulator (Eq. 1), a superior performance is achieved over the model trained with the standard MSE loss function. The interface loss is also considered in the loss function of the GNN assembly model (Eq. 2) to minimize the noncontinuous concentration gap between the assembled geometries, which enables a continuous global concentration prediction. The results indicate that the physics-based loss function serves as an explainable component of the GNN model and teaches the model to utilize the simulation data more efficiently.
Our method, directly learning the transport mechanism from numerical simulation data in mesh format, makes it flexible to design data-driven DL models and further explore the value of simulation data. The model can also be extended to a general GNN framework for learning other PDE models in any complex geometry. Specifically, we can modify the PDE residuals in the simulator loss function to create a physics-guided data-driven DL model for any PDE solver. Our study also has its limitations, which we will address in our future work. In the current model, we only consider different geometries and boundary conditions as input features, while fixing the simulation parameters. In addition, there are still many more neurite networks with different topologies that are not considered in our dataset. To further improve the performance of our GNN model, we will optimize the training dataset by including different simulation parameter settings and different geometries to generalize its application in neurons of broader morphology. We will also study how to combine GNN with physics-informed neural network so that the geometry information of the data can be fully encoded and utilized in the DL framework. Despite these limitations, our GNN model efficiently and accurately predicts the dynamic concentration distribution within complex neurite networks and provides a powerful tool for further study in this field.

IGA-based material transport simulation on complex neuron geometries
In our previous study, we developed an IGA-based solver to perform dynamic material transport simulation in large and complex neurite networks 19 . IGA is an advanced finite element technique that differs from conventional finite element analysis (FEA) for its direct integration of geometrical modeling with numerical solution. With the same smooth spline basis functions 44 utilized as the basis for both geometrical modeling and numerical solution, IGA can accurately represent a wide range of complex geometry with high-order continuity while offering superior performance over FEA in numerical accuracy and robustness 18 . Due to its many performance advantages, IGA has been adopted in a wide variety of research areas, such as linear elasticity 45 , shell analysis [46][47][48] , cardiovascular modeling [49][50][51][52][53] and fluid-structure interaction [54][55][56] , as well as collocation 57,58 . Truncated T-splines 59,60 were developed to facilitate local refinement over unstructured quadrilateral and hexahedral meshes. Blended B-splines 61 and Catmull-Clark subdivision basis functions 62 were investigated to enable improved or even optimal convergence rates for IGA. Our geometric modeling and IGA simulation pipeline 19 is shown in Supplementary Fig. S4. To reconstruct the geometry of different neurite networks for the simulation, we first obtain their morphologies stored in the SWC format from the NeuroMorpho database 40 . With the skeleton and diameter information provided in the SWC files, we adopt the skeleton-sweeping method 49 to generate the hexahedral control mesh of neurite network, and then build trivariate truncated hierarchical B-splines over it. Regarding the governing equations of the simulation, we generalize the 1D motor-assisted transport model 11 to 3D geometry to accurately account for the actual morphology of the neurite. The model is described as a group of "reaction-diffusion-transport" equations, we have where the open set Ω ⊂ R 3 represents the internal space of the single neurite; c 0 , c + and c − are the spatial concentrations of free, incoming (relative to the cell body; retrograde), and outgoing (anterograde) materials, respectively; D is the diffusion coefficient of free materials; u u u + and u u u − are velocities of incoming and outgoing materials, respectively; k ± and k ± are rates of cytoskeletal filament attachment and detachment of incoming and outgoing materials, respectively; and λ ,λ represent the degree of filament attachment at both ends and are also referred to as the "degree of loading" 11 . Regarding the boundary condition, we assume stable concentrations of free and incoming materials at both the incoming end and the outgoing end and set constant values to λ andλ . In this study, we assume the filament system is unipolar that leads to an unidirectional material transport process and ignore c − , u u u − , k − , k − terms in Eq. 3.
To obtain a physically realistic transport process in 3D, we assume that the flow of transport is incompressible and solve the steady-state Navier-Stokes equation to derive a physically realistic velocity field inside the single neurite ∇ · u u u = 0 in Ω, where the open set Ω ⊂ R 3 represents the incompressible fluid domain, u u u is the flow velocity, p is the pressure, f f f is the given body force per unit volume, ν is the kinematic viscosity, and ⊗ denotes the tensor product. Regarding boundary conditions, we impose non-slip condition at the neurite wall and apply a parabolic profile inlet velocity for each point on the circular cross section as where u i is the inlet transport velocity defined in our material transport model, r is the distance from the center of the circular cross section to the point, and R is the radius of the circular cross section. The direction of the velocity is perpendicular to the inlet cross section. We utilize IGA to solve Eq. 4 and get the velocity field u u u. Based on the unidirectional transport assumption, the velocity field is used as u u u + = u u u to solve Eq. 3. As a result, we obtain the material concentration at every time step stored in the truncated hierarchical B-spline of the neurite network. These simulation results are then collected and processed as the training data for this study.

A review of GNNs
GNN is a machine learning technique that was first proposed to generalize the existing neural networks to operate on the graph domain 63 . It has been widely used to perform graph analysis in social networks 30,64 , traffic networks 65 and physics 66,67 . In the following, we explain the basic idea of GNN by using the original GNN framework 63 .
A graph G is a pair (V , E ), where V is the set of nodes and E is the set of edges. Given the node features X V , the edge features X E and outputs O, the GNN is trained to learn the embedding state S to establish the global mapping between all the features and outputs. Since each node is naturally defined by its features and the related nodes in a graph, the nodal embedding state s v and the nodal output o v can be produced locally as o where are the features of node v, the features of its edges, the embedding states and the features of the nodes in the neighborhood of v, respectively. f t is the local transition function and f o is the local output function. Both f t and f o are the parametric functions including the parameters to be trained and shared among all nodes. Let F t (the global transition function) and F o (the global output function) be stacked versions of f t and f o for all nodes in a graph, respectively. Then, we get a compact form of the global mapping: Based on Banach's fixed point theorem 68 , GNN uses the following iterative scheme to compute the embedding states: where S k is the k-th iteration of S. The training of the aforementioned GNN model is straightforward. With the target outputõ v for the supervision, the loss function can be defined as where p is the number of supervised nodes. Then the gradient-descent strategy is used to update and learn the parameters in f t and f o . The original GNN framework has the limitation that they can only handle the nodes embedded in fixed graphs and have to be re-trained whenever the graph structure changes. To address this limitation, GraphSAGE 30 was proposed to generalize the GNN algorithms to learn the nodes embedded in dynamic graphs. The key idea of GraphSAGE is to learn how to aggregate feature information from a node's local neighborhood. When the aggregator weights are learned, the embedding of an unseen node can be generated from its features and neighborhood. To further generalize the concept of GNN, several different GNN frameworks were proposed to integrate different DL models into one framework, such as the message passing neural network (MPNN) 69 , the non-local neural network (NLNN) 70 and the graph network (GN) 35 . In our study, the extensive geometry of neurite networks contributes to different mesh structures for prediction. In particular, different lengths of the pipe lead to different numbers of cross sections along the pipe skeleton. Therefore, we implement GN in our GNN model to ensure the model is suitable for any neuron geometry.

Model evaluation
Two performance metrics were used to evaluate the accuracy of the predicted concentration distributions from the our algorithm: mean absolute error (MAE) and mean relative error (MRE). For each predicted result, the MAE is defined by where N denotes the number of elements in the output, c P i and c G i denote the predicted and ground truth concentration values of the i th node in a given mesh, respectively. For each predicted result, the MRE is defined by where max c G and min c G denote the maximum and minimum nodal concentration values from the ground truth result, respectively.  Test loss Epoch number Figure S3. The test loss vs epoch curve for the GNN assembly model. Figure S4. An overview of the IGA-based neuron material transport simulation pipeline. Given the geometry information (skeleton and diameter) of the neurite network stored in a SWC file, the spline modeling module first generates allhexahedral control mesh and builds volumetric splines over the mesh. Then, the IGA solver takes the volumetric splines as input and computes the dynamic concentration results. Figure S5. The prediction error of NMO_00865 using the model trained with zebrafish neuron dataset. Logarithmic scale is used to highlight the distribution pattern. The regions with higher prediction error are labeled in red circles where the geometry shows high curvature or sharp radius change.