# Hot Isostatic Pressing (HIP) assisted diffusion bonding between CuCr1Zr and
AISI 316L for application to the Super Proton Synchrotron (SPS) internal beam
dump at CERN

S. Pianese<EMAIL_ADDRESS>A. Perillo Marcone  antonio.perillo-
<EMAIL_ADDRESS>F.-X. Nuiry  M. Calviani<EMAIL_ADDRESS>K. Adam
Szczurek  G. Arnau Izquierdo  P. Avigni  S. Bonnin  J. Busom Descarrega  T.
Feniet  K. Kershaw  J. Lendaro  A. Perez Fontenla  CERN, 1211 Geneva 23,
Switzerland  T. Schubert  IFAM, Winterbergstrasse 28, 01277 Dresden, Germany
S. Sgobba  CERN, 1211 Geneva 23, Switzerland  T. Weißgärber  IFAM,
Winterbergstrasse 28, 01277 Dresden, Germany

###### Abstract

The new generation internal beam dump of the Super Proton Synchrotron (SPS) at
CERN will have to dissipate approximately 270 kW of thermal power, deposited
by the primary proton beam. For this purpose, it is essential that the cooling
system features a very efficient heat evacuation. Diffusion bonding assisted
by Hot Isostatic Pressing (HIP) was identified as a promising method of
joining the cooling circuits and the materials of the dump’s core in order to
maximise the heat transfer efficiency. This paper presents the investigation
of HIP assisted diffusion bonding between two CuCr1Zr blanks enclosing SS 316L
tubes and the realisation of a real size prototype of one of the dump’s
cooling plate, as well as the assessments of its cooling performance under the
dump’s most critical operational scenarios. Energy-dispersive X-ray (EDX)
spectroscopy, microstructural analyses, measurements of thermal conductivity
and mechanical strength were performed to characterize the HIP diffusion
bonded interfaces (CuCr1Zr-CuCr1Zr and CuCr1Zr-SS 316L). A test bench allowed
to assess the cooling performance of the real size prototype. At the bonded
interface, the presence of typical diffusional phenomena was observed.
Moreover, measured tensile strength and thermal conductivity were at least
equivalent to the lowest ones of the materials assembled and comparable to its
bulk properties, meaning that a good bonding quality was achieved. Finally,
the real size prototype was successfully tested with an ad-hoc thermal test
bench and with the highest operational thermal power expected in the new
generation SPS internal beam dump. These results demonstrated the possibility
to use HIP as manufacturing technique for the cooling plates of the new
generation SPS internal beam dump, but they also open up the way for further
investigations on its exploitability to improve the cooling performance of any
future high intensity beam intercepting device or in general devices requiring
very efficient heat evacuation systems.

##  I Introduction

The new generation internal beam dump (Target Internal Dump Vertical Graphite,
TIDVG#5 [1] of the Super Proton Synchrotron (SPS) at CERN, built in the
framework of the Large Hadron Collider (LHC) Injectors Upgrade (LIU) Project
[2] has been installed in the long straight section 5 (LSS5) of the SPS during
the CERN’s Long Shutdown 2 (LS2), 2019-2020. Higher beam intensities and
energies, resulting from the LIU upgrades, are making the cooling of this
device, and in general of most beam intercepting devices at CERN, increasingly
demanding. A key problem is obtaining high levels of heat transfer efficiency
between the core beam absorbing materials and the cooling circuits.
Specifically, the SPS internal beam dump is required to evacuate a maximum
average beam thermal power of approximately 270 kW, more than 4 times the
limit of its predecessor TIDVG#4 [3]. As shown in Figure 1, the core of the
SPS internal dump is made of an array of absorbing blocks consisting in 4.4 m
of isostatic graphite, 0.2 m of TZM and $\sim$0.4 m of pure tungsten, enclosed
within 2.5 m long CuCr1Zr cooling plates and under ultra-high vacuum
conditions (UHV). Whenever high energy proton beams need to be dumped (i.e. in
case of emergency, during LHC beam setup or LHC filling, machine developments
and to dispose of the part of the beam for fixed targets remaining after the
slow-extraction process), they are deflected onto the absorbing blocks by an
upstream set of three vertical and three horizontal kickers. The CuCr1Zr
cooling plates must dissipate most of the thermal power carried by the proton
beam and diffused by conduction from the absorbing blocks to the CuCr1Zr heat
sinks. As a consequence, their heat evacuation efficiency is crucial for the
survivability of the absorber system.

Figure 1: Cross section of the new generation SPS internal dump core, TIDVG#5

In the previous generation device [4], the CuCr1Zr plates were mechanically
clamping stainless steel tubes. However, interfaces play a major role in the
heat evacuation mechanism as they represent a resistance to the heat transfer.
For this reason, a possible way to obtain the highest possible thermal contact
conductance was identified by means of diffusion bonding. Hot Isostatic
Pressing assisted diffusion bonding was previously identified as promising
joining technique between copper alloys and stainless steel [5]; by combining
elevated temperatures (up to 2000°C) and isostatic gas pressure (up to 200
MPa), HIP is mainly used to heal casting or additive manufacturing defects and
for powder metallurgy parts consolidations [6, 7, 8]. Although HIP has been
extensively studied as joining technique for several materials [9, 10, 11], it
has not been used so far to manufacture by diffusion bonding large cooling
plates for ultra-high vacuum (UHV) applications. The core of the new
generation SPS beam dump is 5 m long, but, given the maximum dimensions of HIP
ovens available, the HIPed cooling plates were divided in four parts of 2.5 m
each: two bottom upstream/downstream and two top upstream/downstream (Figure
1). The objective of this work was to characterize the HIP diffusion-bonded
interfaces between CuCr1Zr and SS 316L, as well as between CuCr1Zr and
CuCr1Zr, and to assess the cooling performance of a real size prototype of
cooling plate under the dump’s most critical operational scenario.

The use of HIPing to diffusion bond CuCr1Zr with SS316L has also been used
with success for the CERN Proton Synchrotron (PS) Internal dump [12] as well
as for the Beam Dump Facility target prototype beam test [13] to diffusion
bond dissimilar refractory metals. Because of increasingly demanding particle
beams, this manufacturing technique could be of interest to improve the
cooling performances of many other beam intercepting devices (BIDs) such as
the new ISOLDE beam dumps, planned to be refurbished by 2025 within the scope
of the EPIC Project [14], the new 30 kW beam dump of the J-PARC Main Ring (MR)
[15], as well as absorbing devices within the Proton Improvement Plan-II (PIP-
II) Project [16] at FNAL, the International Linear Collider (ILC) photon dumps
or in general any other device requiring very efficient heat evacuation
systems.

##  II Methodology

The needed components for a real size cooling plate prototype were
manufactured, assembled at CERN and diffusion bonded by a HIP cycle. A sub-
part of the prototype was cut and used to manufacture samples to characterize
the diffusion bonded interfaces by means of microstructural analyses, Energy-
dispersive X-ray (EDX) spectroscopy and measurements of thermal conductivity
and mechanical strength. Measurements of electrical conductivity and Brinell
hardness were performed on the bulk of CuCr1Zr to quantify the evolution of
these properties after the HIP cycle and the following thermal treatments. The
rest of the prototype was used within a built-in-house test bench to assess
its cooling performance under the dump’s most critical operational scenario.

###  II.1 Materials

The prototype to be assembled by HIP diffusion bonding (Figure 2) consisted of
two precipitation hardened CuCr1Zr blanks from ZOLLERN GMBH (Germany), 2.7 m
long, embedding in a sandwich two U-bent SS 316L (EN 1.4435) tubes. Table 1
and Table 2 show the chemical composition of this copper alloy and the SS
316L, respectively. The SS 316L cooling tubes are cold bent and annealed. The
annealing cycle for SS 316L consists in a heating-up phase at 200-300 °C/h up
to 950 °C, followed by a plateau at this temperature of at least 2h and
natural cooling-down to ambient temperature. The furnace employed was of the
type “Grand-Four” UGINE for thermal treatments under vacuum, manufactured by
Alstom in 1974. It has an exploitable inner diameter of 990 mm for 6000 mm
height and a maximum capacity of 1000 kg. The maximum reachable temperature is
1100 °C with heating rates ranging from 50 °C/h to 300 °C/h. The cooling is
for natural convection. The CuCr1Zr blanks are multi-directionally forged
according to the standard EN 12420 and heat treated as described below. Refs.
[17, 18] indicate that the following thermal treatments sequence, applied on
CuCr1Zr, allows to obtain the highest values of electrical and thermal
conductivity, hardness, tensile properties, good ductility and machinability:

  * •

Solution annealing between 950-1000 °C for 30 min followed by water quenching;

  * •

Precipitation hardening through ageing at 480-500 °C for 2 h and slow cooling
in air.

Hereinafter in the paper, the material shall be considered to have followed
this treatment and it will be referred as “precipitation hardened CuCr1Zr” or
simply “CuCr1Zr”.

After the HIP cycle, despite the high temperatures, the grains’ structure does
not change, most likely due to the hydrostatic pressure at which the material
is exposed during treatment, but the thermal and mechanical properties need to
be recovered using the same treatment. Table 3 shows how electrical
conductivity and Brinell hardness, measured on the rear side of the prototype,
change after HIP and following the different steps of the heat treatment
cycle. For this application, thermal conductivity was slightly more important
than mechanical strength, therefore the aging phase was performed closer to
the upper limit of 500 °C. Indeed, a higher temperature promotes the
precipitation but induces also a slight coarsening of the precipitates thus
lowering the mechanical strength.

Table 1: CuCr1Zr chemical composition [wt.%] and material designation according to EN 12420 and product analysis Material designation |  |  Element |  Cu |  Cr |  Fe |  Si |  Zr |  Others total
---|---|---|---|---|---|---|---|---
Symbol |  Number |  EN12420 |  Min/Max |  Rem - |  0.5/ 1.2 |  -/ 0.8 |  -/ 0.1 |  0.03/ 0.3 |  -/ 0.2
CuCr1Zr |  EN: CW106C; UNS: C18150 | Analysis product |  Rem |  0.8 |  $<$.01 |  $<$.01 |  .09 |  -
Table 2: SS316L chemical composition [wt.%] and material designation according to EN10088-3 and product analysis Material designation |  |  C |  Si |  Mn |  P |  S |  N |  Cr |  Cu |  Mo |  Ni
---|---|---|---|---|---|---|---|---|---|---|---
AISI 316L |  EN: 10088-3 |  0.030 |  1.00 |  2.00 |  0.045 |  0.030 |  0.11 |  17.0/ 19.0 |  - |  2.50/ 3.00 |  12.5/ 15.0
Product analysis |  0.024 |  0.42 |  1.54 |  0.029 |  0.008 |  0.045 |  17.59 |  - |  2.56 |  13.08
Table 3: Electrical conductivity and Brinell hardness evolution with heat treatments after HIP Heat treatment |  Electrical conductivity |  Hardness
---|---|---
HIP |  41.7 MS/m |  45 HB
Solution annealing and water-quenching |  25 Ms/m |  60 HB
Ageing |  51 MS/m |  108 HB

###  II.2 Prototype design

Figure 2 shows the overall dimensions of the two precipitation hardened
CuCr1Zr halves assembled with the tubes, corresponding to 227 mm (W) × 116 mm
(H) × 2700 mm (L).

Figure 2: Design details of the real size cooling plate prototype.

The tubes inner and outer diameters are 15 mm and 18 mm, respectively. The
rear part of the assembly, containing a straight tube of SS 316L, is meant to
be cut off to be used for the characterization of the bonded interfaces after
the HIP cycle and heat treatments. Before the HIPing cycle, to ease as much as
possible the diffusion bonding at the interfaces, the parts have been cleaned.
The cleaning consists in de-greasing the parts by immersion in a solution of a
compatible alkaline detergent (pH 9.7) at 60°C for 60 minutes, while an
ultrasonic agitation is set for 10 to 20 minutes, and final rinsing with
demineralized water followed by drying with clean compressed air or nitrogen.

###  II.3 HIP capsule

During the HIP cycle, pressure is hydrostatically applied inside the tubes and
outside the precipitation hardened CuCr1Zr blanks. In order to promote
diffusion bonding at the components’ interfaces, these latter must be kept
under vacuum during the HIPing cycle. As shown in Figure 3, a capsule encloses
the assembly and keeps the vacuum leak tightness all over the cycle duration.
The HIP capsule is made of six SS 304L plates: the front and rear plates are 8
mm thick whereas the rest is 2 mm. The plates are welded together, the lateral
and top/bottom by edge welds, while the front and rear by lap welds with the
rest of the plates. In order to lower as much as possible the stresses on the
welds during the HIP cycle, the capsule dimensions are tailored on the inner
CuCr1Zr/tubes assembly. The eight corners formed by the six SS 304L plates are
filled with “plugs”, which are melted during the welding process making more
uniform beads. All the welds of the capsule are tested by dye penetrant and
total helium leak test. The dye penetrant tests are performed according to the
standard ISO 3452-1-2013 and all the welds should be conform at least to the
acceptance level 2X of the standard ISO 23277-2009-EN-BS. The penetrant
liquid, the developer and cleaner used are respectively Ardrox penetrant 907
PB (red dye), Ardrox 9D1 B and Ardrox 9 PR 88. The total helium leak test is
performed according to the standard EN 1779 and NF EN ISO 20485. The capsule
is evacuated and connected to a detector type Leybold L300 or ASM 142
calibrated on helium-4 ( the gas tracer). A higher concentration of the gas
tracer is created around the capsule by enclosing it in a bag. The leak rate
is calculated based on the signal recorded and the helium-4 concentration
inside the bag. The capsule is accepted if there are no leaks within the
sensibility of the test, $1E^{-9}mbar.l/s$.

Figure 3: HIP capsule components and main characteristics.

It is important to reach a good vacuum level (at least $10^{-3}mbar$) in order
to avoid that air pockets remain trapped at the interfaces thus potentially
lowering the bonding quality. Before the HIP cycle, the capsule is pumped down
by means of the lateral tube shown in Figure 4. CuCr1Zr, SS 304L and SS 316L
have very similar thermal expansion coefficients within the operational HIP
temperatures, ranging from $15*10^{-6}m/mK^{-1}$ at room temperature up to
almost $20*10^{-6}m/mK^{-1}$ at 1000 °C: stresses due to relative thermal
expansions and contractions are then limited.

Figure 4: The HIP capsule containing the CuCr1Zr and tubes assembly, PT and
He-leak tested with the lateral tube for pumping

###  II.4 HIPing cycle

The capsule is transported into a HIP unit: a system applying high pressure
with argon gas and high temperature. Table 4 shows the main characteristics of
the HIP unit used.

Table 4: HIP unit characteristics Manufacturer |  Model |  Max load dimensions [mm] |  Max load weight [kg] |  Max Temperature [°C] |  Max pressure [MPa]
---|---|---|---|---|---
ASEA |  QIH307 ”Mega-HIP” |  $\phi$ 1300 h 3200 |  10000 |  1240 |  140

As shown in Figure 5, the HIP cycle phases and the parameters chosen are the
following:

  * •

Phase 1,“$Starting-phase$”: the pressure is raised up to 22.5 MPa and the
temperature is kept constant to “ambient” value (around 50 °C). The parameters
of this phase might change as they depend on the specific HIP unit used.

  * •

Phase 2, “$Heating-phase$”: linear pressure and temperature increase up to 105
MPa and 950 °C. The temperature increases at a rate of $\sim$ 5 °C/min.

  * •

Phase 3, “$Sustaining-phase$”: the following parameters of temperature and
pressure are kept over 180 min:

    * –

Temperature: 950 °C

    * –

Pressure: 105 MPa

  * •

Phase 4, “$Cooling-phase$”: pressure and temperature are decreased to ambient
value. The temperature decreases by natural cooling at a rate of $\sim$ 5
°C/min.

Figure 5: Temperature and pressure over time of the HIP cycle performed for
the prototype of cooling plate of the SPS beam dump

###  II.5 Ultrasonic testing examination and capsule opening

Inspection by Ultrasonic Testing (UT) is a preliminary and fast method to
check if the parts have achieved an intimate contact after the HIP cycle. The
interfaces that can be easily inspected are the SS 304L capsule/CuCr1Zr and
CuCr1Zr/CuCr1Zr as shown in Figure 6. If the HIPing process is successful, the
different parts of the assembly are in intimate contact and there are very few
or no reflections at the interfaces. Specifically, at the interface between
the SS 304L capsule and CuCr1Zr, the ultrasound beam will only be reflected in
a small portion due to the slight difference in acoustic impedance of the
joined materials. This interface was inspected with an ALPHA2, 15 MHz and 6 mm
diameter probe. With a gain of 40 dB and the speed of sound in stainless steel
of 5800 m/s, the reflection obtained was less than 5%. The interface between
CuCr1Zr/CuCr1Zr, given its bigger depth, was inspected with the lower
frequency SEB4 probe of 4 MHz and 20 mm diameter. The only reflection obtained
was the back-wall echo corresponding to the thickness of the assembly. This
inspection was carried out by setting a gain of 60 dB and the speed of sound
in copper of about 4800 m/s. The ultrasound beams generator was an USM Go+.
The probes and the UT generator were provided by General Electric. After the
UT inspections, the capsule was removed by machining.

Figure 6: Interfaces inspected by ultrasounds after HIP cycle

###  II.6 HIP diffusion bonded interface examination

The interfaces were inspected with optical and electronic microscopes to
assess the diffusion bonding quality. The equipment used for the examinations
consisted in:

  * •

Digital microscope KEYENCE VHX 6000;

  * •

Scanning Electron Microscope (SEM), field emission gun FEG Sigma (ZEISS) with
InLens (Secondary Electron), Everhart-Thornley Secondary Electron (SE2), and
back-scattered electron (AsB) detectors for imaging;

  * •

$50mm^{2}$ X Max Electron dispersive X-ray spectroscopy (EDS) detector
(Oxford), AZTEC software.

Figure 7 shows that the interface CuCr1Zr - SS 316 L, is free from major
defects. The bonding zone shows a low amount of micro-porosity on the CuCr1Zr
side, most probably corresponding to the Kirkendall effect [19] due to the
different diffusion coefficients between the bonding materials, a first
indicator that the diffusion took place.

Figure 7: SEM image at 3000 X and 15kV accelerating voltage. A continuous
$\delta$-ferrite layer is visible on the SS side, whereas micro porosities due
to the Kirkendall effects and Zr inclusions are distinguishable on the CuCr1Zr
side. Figure 8: EDX elemental line-scans acquired at 20 keV accelerating
voltage in a direction perpendicular to the interface in order to assess the
diffusion phenomena. The “x” axis represents in $\mu$m the yellow line at the
bottom of the picture and the “y” axis gives the concentration of the chemical
elements along this line. On the SS 316L side, there is a clear excess of Cr
up to 30 Wt% and a lack of Ni with some diffused Cu. On the CuCr1Zr side, it
can be distinguished a peak corresponding to a Zr inclusion and some diffused
Fe.

As shown by Figure 8, several EDX analyses of the interface CuCr1Zr – SS 316L
revealed a general trend: a few micrometres Cr rich and Ni poor layer with
some diffused Cu is present on the steel side whereas Zr rich precipitates
with some diffused Fe are on the CuCr1Zr side [20]. The darker continuous
layer of about 5 µm thickness visible on the steel side is $\delta$-ferrite.
According to the Schaeffler diagram for stainless steel, nickel loss at the SS
316L/CuCr1Zr interface, Cr diffusion to SS 316L, plus Cr enrichment up to 30
wt% are responsible for the formation of a fully ferritic layer at the joint.
Figure 9 also shows large dark Cr and Zr precipitates consolidated during the
aging process which diminish in a layer of about 50 µm close to the SS316L
/CuCr1Zr interface.

Figure 9: SEM image 500X of the CuCr1Zr-SS 316L interface. Dark “Cr” and “Zr”
precipitates consolidated during aging are visible on the CuCr1Zr side

Figure 10 shows an optical microscopy of the interface CuCr1Zr-CuCr1Zr, which
appears clearly discernible, homogeneous and without significant voids
population. Interfacial grain boundaries lie in the plane of the joint.

Figure 10: Optical microscopy of a sample containing the HIP diffusion bonded
CuCr1Zr-CuCr1Zr interface. The surface inspected has been mechanically
polished and chemically etched.

###  II.7 Thermal conductivity measurements

Figure 11 shows the measurements of thermal conductivity, from room
temperature to 300 °C, of pure copper, precipitation hardened CuCr1Zr, SS 316L
and on samples containing their respective HIP diffusion-bonded interfaces,
performed by using a NETZSCH nanoflash LFA447.

Figure 11: Thermal conductivity measurements at the materials’ bulks and
interfaces.

At the interface CuCr1Zr-CuCr1Zr (dark blue line), the thermal conductivity is
as high as that measured in the bulks of precipitation hardened CuCr1Zr
(orange line) and pure copper (light blue line), with values ranging between
$350\pm 10Wm^{-1}K^{-1}$ and $360\pm 10Wm^{-1}K^{-1}$. Similarly, the thermal
conductivity at the interface CuCr1Zr-SS 316L (green line) and the bulk of SS
316L (yellow line) are close and fall in the range between $38\pm
1Wm^{-1}K^{-1}$ and $45.0\pm 1.4Wm^{-1}K^{-1}$. These results indicated that
at the HIP diffusion bonded interfaces the values are limited by the material
with the lowest thermal conductivity in the joint, however higher than the
latter (see e.g. CuCr1Zr-SS 316L interface).

###  II.8 Mechanical strength measurements

Tensile tests at room temperature, according to the DIN EN ISO 6892-1, of the
HIP diffusion bonded interfaces were performed in order to quantify their
mechanical strengths. Samples made of CuCr1Zr and SS 316L as well as
containing the HIP diffusion bonded interfaces were prepared by electro-
erosion. The initial cross section dimensions and gauge length are the
following (Figure 12):

  * •

Width: $b_{0}=1.5mm$

  * •

Thickness $a_{0}=2mm$

  * •

Cross-sectional area: $S_{0}=3mm^{2}$

  * •

Gauge length: $L_{0}=4.6mm$

Figure 12: Tensile samples dimensions

Figure 13 shows the tested tensile samples containing the CuCr1Zr-CuCr1Zr and
CuCr1Zr-SS 316L interfaces within the gauge length.

Figure 13: Tensile samples containing (a) the CuCr1Zr-CuCr1Zr interface and
(b) the CuCr1Zr-SS 316L interface.

For the CuCr1Zr-SS 316L samples, the fracture appears to be ductile and it is
not at the interface but within CuCr1Zr. The measured mean tensile strength of
360 $\pm$ 2 MPa is comparable with precipitation hardened CuCr1Zr, as it is
shown by the dotted line in Figure 14. Similarly, Figure 15 shows the results
of the tensile samples containing the CuCr1Zr-CuCr1Zr interface. The mean
measured tensile strength of 348 $\pm$ 2 MPa is very close to the bulk
precipitation hardened CuCr1Zr of 352 $\pm$ 3 MPa. The presence of a diffusion
bonded interface does not seem to affect the mechanical strength, which is,
instead, limited by the lowest one of the joining materials. Nevertheless, the
samples containing the CuCr1Zr/SS 316L interface exhibit a lower elongation at
necking with respect to the bulky CuCr1Zr ones. This difference is smaller for
the CuCr1Zr/CuCr1Zr samples. However, the total elongations are comparable and
in most cases close to 30%. Table 5 summarizes the mean values of tensile
strengths.

Table 5: Mean tensile strengths of the specimens with the interfaces CuCr1Zr/CuCr1Zr, CuCr1Zr/SS 316L and for specimens of CuCr1Zr without interface. Samples | Tensile strength [MPa]
---|---
CuCr1Zr/SS 316L | $360\pm 2$
CuCr1Zr/CuCr1Zr | $348\pm 2$
Bulk CuCr1Zr | $352\pm 3$
Figure 14: Room temperature tensile stress-strain curves for CuCr1Zr-SS 316L
specimens and bulk hardened CuC1Zr as a reference. Figure 15: Room temperature
tensile stress-strain curves for CuCr1Zr-CuCr1Zr specimens and bulk hardened
CuCr1Zr.

###  II.9 Cooling performance test bench

In order to assess the ability of the HIP manufactured prototype to evacuate
the thermal power deposited in the most critical regions of the SPS internal
dump’s cooling plates and in the worst operational scenario, as well as to
verify that the cooling performance all over the prototype was homogeneous, a
dedicated test bench was designed. In order to validate the first point, the
cooling plate prototype needed to be able to reach a steady state with an
average temperature well below 300 °C. Indeed, as Figure 16 shows, starting
from 300 °C the measured yield strength of precipitation hardened CuCr1Zr
becomes too low for a robust and reliable operation of the SPS internal dump.

Figure 16: Yield and ultimate strength of precipitation hardened CuCr1Zr from
room temperature to 400°C. Every CuCr1Zr component of the SPS internal dump
shall operate at a temperature always below the limit of 300°C.

As for the homogeneity of the cooling performance, the amount of power
evacuated by the two tubes as well as the cooling rates measured at different
locations of the prototype plate shall be the same within the measurement
errors.

####  II.9.1 Experimental setup

The test bench is made of four main systems/part: the heating and cooling
systems, the monitoring system and the part tested. The heating system
reproduces the thermal power deposited by the particle beam. The cooling
system extracts this power in the same way the actual system would operate on
the dump’s cooling plates in operation and finally the monitoring system gives
feedback about the temperatures and the inlet/outlet flow rates of the part
being tested. Figure 17 shows the main test bench components used to assess
the cooling performance of the HIP manufactured cooling plate prototype.
Specifically, the thermal power deposited by the particle beam is reproduced
by six oxygen-free high thermal conductivity (OFHC) copper heating blocks.
Every heating block has 7 electrical heating cartridges 2 kW each, for a total
power of 84 kW.

Figure 17: Components of the cooling test bench and the respective
configuration.

In order to improve the heat flow, a copper-based high temperature lubricant
with very good thermal conductivity (ORAPI TM900S) couples the heating
cartridges with the heating blocks and these latter with the cooling plate
prototype. The heating blocks are pressed against the prototype plate
establishing an average calculated contact pressure of about 63 kPa. The
temperature distribution at the interface heating blocks - cooling plate is
monitored by 36 PT100s sensors. The whole assembly is enclosed within
insulating fiberglass plates 50 mm thick.

Figure 18: P&ID cooling performance test bench

The prototype plate is equipped with two cooling circuits shown in Figure 17,
following the same design of the SPS internal dump cooling loop. Figure 18
shows the P&ID of the circuits. An inlet flow rate of 5 $m^{3}/h$ and 7 $bar$
of demineralized water is equally split in the two circuits of the prototype.
At the inlet, the temperature is kept constant at 27 °C and it is monitored by
a temperature indicator and transmitter installed directly in the line.
Flowrates and temperatures at the outlet circuits are measured by two
flowmeters type FC 100 LQ from Flow Vision. The electrical heating cartridges,
the PT100 temperature sensors, the flowmeters and the inlet temperature
transmitter are all connected to a PLC and controlled via dedicated software.

####  II.9.2 Test procedure

Every test was carried out by the following steps:

  * •

The cooling circuits are activated;

  * •

The thermal power is gradually increased until reaching the testing level;

  * •

The test bench is left to reach the steady state where the average cooling
plate temperature and the cooling circuit’ s inlet and outlets temperatures
are recorded;

  * •

The power is shut down and the cooling rates at different points of the plates
are measured.

####  II.9.3 Results

The highest tested thermal power was 84 kW, which is roughly equivalent to the
maximum expected average power deposited in the most critical of the four
plates composing the TIDVG#5. Figure 19 shows the increase of average
temperature of all the PT100s sensors over time, from the beginning up to the
steady state. Given the cooling system in place and the amount of power
tested, it took about 32 minutes to reach the steady state at an average
temperature of 136 °C. Figure 20 shows the temperature over time at the inlet
and outlets of the cooling circuits. At steady state condition, this
information along with the flow rates allow to calculate the thermal power
extracted by the circuits. The outlets’ temperatures difference less than 2 °C
suggests that the two circuits behave in the same way and evacuate a similar
amount of thermal power, precisely, 39 kW and 35 kW summing up to 74 kW.
Around 10 kW are lost due to the not perfect insulation of the test bench.
When the power is shut down and the cooling plate temperature starts
decreasing, the cooling rates in several points of the prototype were
measured. Figure 21 shows that the maximum difference between the cooling
rates is below 0.5 °C/s. This result, along with the similar thermal power
evacuated, suggests that the cooling performance of the prototype are
homogeneous.

Figure 19: Average PT100s temperature over time to steady state of the heating
phase. Figure 20: Inlet and outlets water temperature versus time during
heating phase at $84kW$. Figure 21: Maximum cooling rates measured in some
points of the prototype’s surface.

##  III Discussion

The objective of the present study was to carry out a complete
characterization of the CuCr1Zr-CuCr1Zr and CuCr1Zr-SS 316L HIP diffusion
bonded interfaces and to assess the cooling performance of the cooling plate
prototype under the worst operational scenario of the SPS beam internal dump.
In its most narrow definition, diffusion bonding is a process that produces
solid-state coalescence between two clean surfaces. The process is composed of
a three-stage metallurgical phases citeASM:

  1. 1.

The contact area grows to a large fraction of the joint area by localized
deformation of the contacting surface asperities. At the completion of this
stage, the interface boundary is no longer a planar interface, but consists of
voids separated by areas of intimate contact. In these areas of contact, the
joint becomes equivalent to a grain boundary between the grains on each
surface.

  2. 2.

All of the voids in the joints shrink and most are eliminated. The interfacial
grain boundary migrates out of the plane of the joint to a lower-energy
equilibrium. Bonding could be considered essentially complete after this
stage.

  3. 3.

Diffusional processes caused the shrinkage and elimination of voids, now the
only possible path is through the volume of the grains themselves.

For the CuCr1Zr-SS 316L interface, there was a good contact of the surfaces
and two diffusional phenomena also observed in literature [5]:

  * •

Low amount of micro-porosity on the CuCr1Zr side, most probably due to the
Kirkendall effect;

  * •

The Cr rich and Ni poor diffusion layer with some diffused Cu on the SS 316L
side and Cr and Zr rich inclusions with some diffused Fe on the CuCr1Zr side.

The CuCr1Zr-CuCr1Zr interface was clearly discernible, homogeneous and without
significant voids population. For both interfaces, the interfacial grain
boundaries remained in the plane of the joint, therefore the diffusion bonding
cannot be considered fully complete according to the above definition.
Nevertheless, the measurements of thermal conductivity and mechanical
strength, for the CuCr1Zr-SS 316L interface, showed values limited by the SS
316L and comparable to its bulk whereas for the CuCr1Zr-CuCr1Zr interface the
values were comparable to the bulk of precipitation hardened CuCr1Zr. For
these reasons, it can be concluded that the diffusion bonding reached advanced
stage (phase 2). The data obtained from the test bench showed that the
manufactured real size prototype is well capable of withstanding the thermal
power that the most critical of the SPS dump’s cooling plates would undergo in
the worst operational scenario. Finally, the homogeneity of the cooling
performance suggests that the HIP diffusion bonding process, for the
corresponding interfaces all over the prototype plate, is homogeneous.

##  IV Conclusions

The future CERN SPS internal dump, TIDVG#5, will have to cope with higher beam
intensities and repetition dumping rates starting from 2021. The maximum
average thermal power carried by the beam can be as high as approximately 270
kW and about 90% of it is absorbed by the dump and shall be evacuated by the
cooling plates. Their cooling performance is, therefore, crucial. Diffusion
bonding of stainless steel 316L and precipitation hardened CuCr1Zr by HIPing
was thoroughly investigated as manufacturing technique for the dump’s cooling
plates to reduce interfacial thermal resistance between the different
components. Micro-structural analyses and measurements of thermal conductivity
and mechanical strength suggest that diffusion bonding at the interfaces
CuCr1Zr-SS316L and CuCr1Zr-CuCr1Zr reached an advanced stage. Moreover, a test
bench allowed to measure the cooling performance of the real size prototype
with the worst operational scenario and proved that the cooling performance
are uniform and the manufacturing technique can be used for the real cooling
plates of the future SPS internal dump. In conclusion, the present work opens
up the way for further investigations on the exploitability of this technique
for other challenging cooling applications in the field of high energy
physics.

## References

  * Pianese  _et al._ [2018] S. Pianese _et al._ , Design of the future high energy beam dump for the CERN SPS, in  _Proc. 9th Int. Particle Accelerator Conf. (IPAC’18)_ , Vancouver, Canada (2018) pp. 2612–2615, paper WEPMG004.
  * Coupard  _et al._ [2014] J. Coupard, H. Damerau, A. Funken, R. Garoby, S. Gilardoni, B. Goddard, K. Hanke, A. Lombardi, D. Manglunki, M. Meddahi, B. Mikulec, G. Rumolo, E. Shaposhnikova, and M. Vretenar,  _LHC Injector Upgrade_ , Technical Design Report, Vol 1: Protons (Geneva, 2014) CERN-ACC-2014-0337.
  * Rios-Rodriguez  _et al._ [2017] P. Rios-Rodriguez _et al._ , Analysis and operational feedback of the new design of high energy beam dump in the CERN SPS, in  _Proc. 8th Int. Particle Accelerator Conf. (IPAC’17)_ , Copenhagen, Denmark (2017) pp. 3524–3527, paper WEPVA110.
  * Ross [2000] M. Ross,  _The millenium dump TIDVG commissioning - april 2000_, Tech. Rep. SL-Note-2000-050-BT (CERN, Geneva, 2000).
  * Marois  _et al._ [1996] G. L. Marois _et al._ , Hip’ing of copper alloys to stainless steel, Journal of Nuclear Material 233-237 , 927 (1996).
  * Atkinson and Davies [2000] H. Atkinson and S. Davies, Fundamental aspects of hot isostatic pressing: An overview, Metallurgical and material transactions Volume 31 A , 2981 (2000).
  * Kumar  _et al._ [2018] A. Y. Kumar _et al._ , The effects of hot isostatic pressing on parts fabricated by binder jetting additive manufacturing, Additive Manufacturing, Volume 24 , 115 (2018).
  * Wu  _et al._ [2017] J. Wu _et al._ , Effect of hot isostatic pressing loading route on microstructure and mechanical properties of powder metallurgy ti2alnb alloys, Journal of Materials Science & Technology, Volume 33, Issue 2 , 172 (2017).
  * [9] J. B. Descarrega, M. Calviani, T. Hutsch, _et al._ , Application of hot isostatic pressing (HIP) technology to diffusion bond refractory metals for proton beam targets and absorbers at CERN, Mat Design Process Comm. 2020;2:e101. https://doi.org/10.1002/mdp2.101 .
  * [10] N. Ordás _et al._ , Mechanical and microstructural characterization of HIP joints of a simplified prototype of the ITER NHF first wall panel, http://dx.doi.org/10.1016/j.fusengdes.2017.05.087 0920-3796/© 2017 Elsevier B.V .
  * [11] S. Makimura _et al._ , Present status of construction for the muon target in j-parc, Nuclear Instruments and Methods in Physics research A 600 (2009) 146-149, Elsevier B.V, doi:10.1016/j.nima.2008.11.110 .
  * Romagnoli  _et al._ [2018] G. Romagnoli _et al._ , Engineering design and prototyping of the new LIU PS internal beam dumps, in  _Proc. 9th Int. Particle Accelerator Conf. (IPAC’18)_ , Vancouver, Canada (2018) pp. 2600–2603, wEPMG001.
  * [13] E. L. Sola, M. Calviani, _et al._ , Beam impact tests of a prototype target for the beam dump facility at CERN: Experimental set-up and preliminary analysis of the online results, Physical review accelerators and beams 22, 2019, pp. 123001-1 – 123001-16, https://doi.org/10.1103/PhysRevAccelBeams.22.123001 .
  * Catherall  _et al._ [2019] R. Catherall, T.Giles, and G. Neyens, Exploiting the potential of ISOLDE at CERN (the EPIC project), in  _Proc. 10th Int. Particle Accelerator Conf. (IPAC’19)_ , Melbourne, Australia (2019) pp. 3706–3708, paper THPGW053, http://accelconf.web.cern.ch/ipac2019/papers/thpgw053.pdf.
  * Shirakata  _et al._ [2019] M. J. Shirakata, H. Kuboki, and J. Takano, Radiation design of new 30 kw beam dump of J-PARC main ring, in  _Proc. 10th Int. Particle Accelerator Conf. (IPAC’19)_ , Melbourne, Australia (2019) pp. 2005–2008, paper TUPTS032, http://accelconf.web.cern.ch/ipac2019/papers/tupts032.pdf.
  * Lebedev [2015] V. Lebedev,  _The PIP-II Reference Design Report_ , Tech. Rep. (Fermi National Accelerator Laboratory, 2015) https://lss.fnal.gov/archive/design/fermilab-design-2015-01.pdf.
  * Ivanov  _et al._ [2002] A. D. Ivanov _et al._ , Effects of heat treatments on the properties of CUCR1ZR alloys, Journal of Nuclear Mareials, 307-311 , 673 (2002).
  * Holzwarth  _et al._ [2000] U. Holzwarth _et al._ , The recovery of tensile properties of CUCR1ZR alloy after hot isostatic pressing, Fusion Engineering and Design 51-52 , 111 (2000).
  * Pieraggi [2008] B. Pieraggi, in  _Developments in High Temperature Corrosion and Protection of Materials_ , Metal and Surface Engineering (Woodhead Publishing, 2008) Chap. 2, pp. 9–35, diffusion and solid state reactions.
  * Fagnoni  _et al._ [2018] F. Fagnoni, J. B. Descarrega, and E. R. Castro, Characterization of diffusion bonded cucr1zr-cucr1zr and cucr1zr-steel interfaces for PS dump prototype (2018), EDMS document: 2136203, CERN, Geneva.


# Deep learning of material transport in complex neurite networks

Angran Li  Department of Mechanical Engineering, Carnegie Mellon University,
Pittsburgh, PA, USA  Amir Barati Farimani  Department of Mechanical
Engineering, Carnegie Mellon University, Pittsburgh, PA, USA  Department of
Biomedical Engineering, Carnegie Mellon University, Pittsburgh, PA, USA
Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh,
PA, USA  Department of Machine Learning, Carnegie Mellon University,
Pittsburgh, PA, USA  Yongjie Jessica Zhang  Department of Mechanical
Engineering, Carnegie Mellon University, Pittsburgh, PA, USA  Department of
Biomedical Engineering, Carnegie Mellon University, Pittsburgh, PA, USA
<EMAIL_ADDRESS>
###### Abstract

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\sim 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\sim 330$ times faster
compared to IGA simulations.

## Results

### 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_{S}^{p}$) and the bifurcation ($G_{S}^{b}$), 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 $x_{mid}$, the model can assemble different
structures and output the final prediction $x_{o}$ on the graph. In sum, our
GNN model consists of

  * •

The GNN simulators for local prediction in pipe and bifurcation structures;
and

  * •

The GNN assembly model for global prediction in the entire neurite network.

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
$\phi_{e}$, $\phi_{v}$ to update the edge and node features, respectively.
$\phi_{e}$ computes the concentration gradient along the edge and the length
of the edge. $\phi_{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.

Input: Graph, $\mathcal{G}=(\mathcal{V},\mathcal{E},x_{v})$

for  _each edge $e_{ij}\in\mathcal{E}$_ do

Compute edge attributes: $x_{e_{ij}}=\phi_{e}(x_{v_{i}},x_{v_{j}})$;

end for

for  _each node $v_{i}\in\mathcal{V}$_ do

Aggregate edge attributes: $x_{e,v_{i}}=\sum_{j}x_{e_{ij}}$;

Compute output: $x^{\prime}_{v_{i}}=\phi_{v}(x_{v_{i}},x_{e,v_{i}})$;

end for

Output: Graph, $\mathcal{G}=(\mathcal{V},\mathcal{E},x^{\prime}_{v})$

Algorithm 1 GN Block

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.

Figure 2: Results of GNN simulators. (A) The architecture of the GNN
simulators adopts a recurrent “GN Block + MLP Decoder" scheme on the input
graph to compute nodal concentration values at each time step from input nodal
features. The input $n_{v}$-node graph is stored in a $n_{v}\times n_{v}$
adjacent matrix. The nodal features include coordinate vector $x$, simulation
parameter vector $p$ and initial concentration values $c^{t_{k}}$ at time step
$t_{k}$. The $L$-step “GN Block" follows Algorithm 1 to compute the
interaction among nodes and generates a series of updated latent graphs with
hidden nodal embeddings, $h^{1}$, …, $h^{L}$. The “MLP Decoder" outputs nodal
concentration values $c^{t_{k+1}}$ at the next time step from the nodal
embedding $h^{L}$ of the final latent graph. (B, C) The concentration
distribution comparison of pipe and bifurcation simulators, respectively. For
each simulator, we plot the prediction and nodal error results at three
different time steps. The average errors are 3.45% and 4.97% for the pipe and
bifurcation simulators, respectively.

To encourage the model to learn the physics underlying the simulation data, we
add the residuals of the governing equation (Eq. 3 in Methods) to the mean
squared error (MSE) loss function as

$\begin{split}\mathcal{L}_{simulator}&=\frac{1}{N}\sum\limits^{N}_{i=1}\left[(c^{P}_{0,i}-c^{G}_{0,i})^{2}+(c^{P}_{\pm,i}-c^{G}_{\pm,i})^{2}\right]+\frac{1}{N}\sum\limits^{N}_{i=1}\left[\dfrac{\partial
c_{0,i}}{\partial
t}-D\nabla^{2}c_{0,i}+(k_{+}+k_{-})c_{0,i}-k^{\prime}_{+}c_{+,i}-k^{\prime}_{-}c_{-,i}\right]^{2}\\\
&+\frac{1}{N}\sum\limits^{N}_{i=1}\left[\dfrac{\partial c_{\pm,i}}{\partial
t}+\bm{u}_{\pm}\cdot\nabla
c_{\pm,i}-k_{\pm}c_{0,i}+k^{\prime}_{\pm}c_{\pm,i}\right]^{2},\end{split}$ (1)

where $c_{0}$, $c_{+}$ and $c_{-}$ are the spatial concentrations; $D$ is the
diffusion coefficient of materials; $\bm{u}_{+}$ and $\bm{u}_{-}$ are
velocities of materials; $k_{\pm}$ and $k^{\prime}_{\pm}$ 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\;\mu m^{2}/s$, $k_{\pm}=1.0\;s^{-1}$,
$k^{\prime}_{\pm}=0.5\;s^{-1}$ and $u_{i}=0.1\;\mu 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 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).

Figure 3: Results of the GNN assembly model. (A) The GNN assembly model takes
the intermediate prediction $X$ from individual GNN simulators as input and
outputs the final concentration prediction $X^{\prime}$. The GNN assembly
components $G_{A}^{p-p}$, $G_{A}^{p-b}$ and $G_{A}^{b-b}$ adopt the “message-
passing" scheme to aggregate the intermediate prediction from the neighbouring
nodes to update the concentration prediction. The green arrows show the
message-passing process during the assembly. (B, C) The predicted
concentration and the nodal prediction error of NMO_66731 at steady state
($t=25\;s$). (D, E) The predicted concentration and the nodal prediction error
of NMO_66748 at steady state ($t=31\;s$). In (B) and (D), the black arrows
point to the inlet of material. Scalar bar: 20 $\mu m$. Unit for color bars:
$mol/\mu m^{3}$.

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 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.

Table 1: MRE comparison between simulators using different loss functions Loss function | MSE | MSE
---|---|---
|  | \+ PDE residuals
Pipe simulator | 10.9 $\pm$ 3.45 % | 6.10 $\pm$ 1.25 %
Bifurcation simulator | 13.5 $\pm$ 2.47 % | 7.20 $\pm$ 0.78 %

### 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_{A}^{p-p}$; (b)
pipe and bifurcation $G_{A}^{p-b}$; and (c) bifurcation and bifurcation
$G_{A}^{b-b}$. 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
$G_{A}^{p-p}$, $G_{A}^{p-b}$ and $G_{A}^{b-b}$ 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

$\mathcal{L}_{assembly}=\frac{1}{N}\sum\limits^{N}_{i=1}\left\\{(c^{P}_{0,i}-c^{G}_{0,i})^{2}+(c^{P}_{\pm,i}-c^{G}_{\pm,i})^{2}\right\\}+\alpha\frac{1}{M}\sum\limits^{M}_{j=1}\left\\{(c^{s_{1},interface}_{0,j}-c^{s_{2},interface}_{0,j})^{2}\right.\left.+(c^{s_{1},interface}_{\pm,j}-c^{s_{2},interface}_{\pm,j})^{2}\right\\},$
(2)

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. $\alpha$ 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.

Figure 4: The concentration prediction of material transport in the neurite
network geometry of NMO_06846, NMO_06840, NMO_112145, NMO_32235 and NMO_32280.
(A-C, G, H) The predicted concentration results of steady state at $t=15\;s$,
$26\;s$, $19\;s$, $28\;s$ and $32\;s$, respectively. The white arrows point to
the inlet of material. Scalar bar: 20 $\mu m$. Unit for color bars: $mol/\mu
m^{3}$. (D-F, I, J) The nodal errors between the ground truth and predicted
concentration. Logarithmic scale is used to highlight the distribution
pattern. Figure 5: The concentration prediction of material transport in the
neurite network geometry of NMO_54504, NMO_54499 and NMO_00865. (A, C, E) The
predicted concentration results of steady state at $t=12\;s$, $41\;s$, and
$63\;s$, respectively. The black arrows point to the inlet of material. Scalar
bar: 20 $\mu m$. Unit for color bars: $mol/\mu m^{3}$. (B, D, F) The nodal
errors between the ground truth and predicted concentration. Logarithmic scale
is used to highlight the distribution pattern.

### 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 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.

Table 2: Statistics of all tested complex neurite networks

Species | Model name | Mesh | Bifurcation | IGA computation | GNN prediction | Speedup ratio | GNN prediction
---|---|---|---|---|---|---|---
(vertices, elements) | number | (nodes, time (mins)) | time (mins) | (IGA vs GNN) | MRE
Zebrafish | NMO_66731 (Fig. 3B) | (127,221, 112,500) | 15 | (8, 468) | 1.6 | 293 | 6.7%
NMO_66748 (Fig. 3D) | (282,150, 249,660) | 35 | (10, 672) | 3.9 | 172 | 7.3%
NMO_06846 (Fig. 4A) | (116,943, 101,880) | 20 | (8, 413) | 2.1 | 197 | 7.4%
NMO_06840 (Fig. 4B) | (280,434, 248,040) | 35 | (10, 705) | 4.4 | 160 | 7.6%
Mouse | NMO_112145 (Fig. 4C) | (110,985, 98,460) | 9 | (8, 350) | 1.3 | 269 | 8.2%
NMO_32235 (Fig. 4H) | (96,714, 85,680) | 9 | (6, 320) | 1.1 | 291 | 8.3%
NMO_32280 (Fig. 4I) | (131,967, 117,000) | 12 | (8, 493) | 1.5 | 329 | 8.1%
NMO_54504 (Fig. 5A) | (116,943, 101,880) | 32 | (8, 436) | 3.3 | 132 | 8.3%
NMO_54499 (Fig. 5C) | (524,871 459,360) | 127 | (20, 759) | 6.1 | 124 | 8.9%
NMO_00865 (Fig. 5E) | (1,350,864, 1,179,900) | 356 | (40, 908) | 7.2 | 126 | 9.2%

  * •

Note: Each node in Bridges has 28 cores.

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 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.

## Methods

### 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

$\left\\{\begin{array}[]{lr}\dfrac{\partial c_{0}}{\partial
t}-D\nabla^{2}c_{0}=-(k_{+}+k_{-})c_{0}+k^{\prime}_{+}c_{+}+k^{\prime}_{-}c_{-}&\text{in
$\Omega$},\vspace{1.5mm}\\\ \dfrac{\partial c_{\pm}}{\partial
t}+\bm{u}_{\pm}\cdot\nabla
c_{\pm}=k_{\pm}c_{0}-k^{\prime}_{\pm}c_{\pm}&\text{in $\Omega$},\\\
c_{0}=c,\indent c_{+}=\lambda c&\text{at incoming end},\\\
c_{0}=\tilde{c},\indent c_{-}=\tilde{\lambda}\tilde{c}&\text{at outgoing
end},\end{array}\right.$ (3)

where the open set $\Omega\subset\mathbb{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; $\bm{u}_{+}$ and $\bm{u}_{-}$ are velocities of
incoming and outgoing materials, respectively; $k_{\pm}$ and
$k^{\prime}_{\pm}$ are rates of cytoskeletal filament attachment and
detachment of incoming and outgoing materials, respectively; and $\lambda$,
$\tilde{\lambda}$ 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
$\lambda$ and $\tilde{\lambda}$. In this study, we assume the filament system
is unipolar that leads to an unidirectional material transport process and
ignore $c_{-}$, $\bm{u}_{-}$, $k_{-}$, $k^{\prime}_{-}$ 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

$\begin{cases}\nabla\cdot\bm{u}=0&\text{in $\Omega$},\\\
\nabla\cdot(\bm{u}\otimes\bm{u})+\nabla p=\nu\Delta\bm{u}+\bm{f}&\text{in
$\Omega$},\end{cases}$ (4)

where the open set $\Omega\subset\mathbb{R}^{3}$ represents the incompressible
fluid domain, $\bm{u}$ is the flow velocity, $p$ is the pressure, $\bm{f}$ is
the given body force per unit volume, $\nu$ is the kinematic viscosity, and
$\otimes$ 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

$u(r)=u_{i}(1-(r/R)^{2}),$ (5)

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 $\bm{u}$. Based on the
unidirectional transport assumption, the velocity field is used as
$\bm{u}_{+}=\bm{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 $\mathcal{G}$ is a pair $(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}$ is the set of nodes and $\mathcal{E}$ is the set of edges. Given
the node features $\textbf{X}_{\mathcal{V}}$, the edge features
$\textbf{X}_{\mathcal{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 $\textbf{s}_{v}$ and the nodal
output $\textbf{o}_{v}$ can be produced locally as

$\displaystyle\textbf{s}_{v}$
$\displaystyle=f_{t}(\textbf{x}_{v},\textbf{x}_{\epsilon},\textbf{s}_{ne[v]},\textbf{x}_{ne[v]}),$
(6) $\displaystyle\textbf{o}_{v}$
$\displaystyle=f_{o}(\textbf{s}_{v},\textbf{x}_{v}),$ (7)

where $\textbf{x}_{v},\textbf{x}_{e},\textbf{s}_{ne[v]},\textbf{x}_{ne[v]}$
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:

S
$\displaystyle=F_{t}(\textbf{S},\textbf{X}_{\mathcal{V}},\textbf{X}_{\mathcal{E}}),$
(8) O $\displaystyle=F_{o}(\textbf{S},\textbf{X}_{\mathcal{V}}).$ (9)

Based on Banach’s fixed point theorem [68], GNN uses the following iterative
scheme to compute the embedding states:

$\textbf{S}^{k+1}=F_{t}(\textbf{S}^{k},\textbf{X}_{\mathcal{V}},\textbf{X}_{\mathcal{E}}),$
(10)

where $\textbf{S}^{k}$ is the $k$-th iteration of S. The training of the
aforementioned GNN model is straightforward. With the target output
$\tilde{\textbf{o}}_{v}$ for the supervision, the loss function can be defined
as

$loss=\sum_{i=1}^{p}||\tilde{\textbf{o}}_{v,i}-\textbf{o}_{v,i}||,$ (11)

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

$MAE=\sqrt{\sum\limits^{N}_{i=1}\frac{1}{N}(c^{P}_{i}-c^{G}_{i})^{2}},$ (12)

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

$MRE=\frac{MAE}{\max\left|\left|c^{G}\right|\right|-\min\left|\left|c^{G}\right|\right|}\times
100\%,$ (13)

where $\max\left|\left|c^{G}\right|\right|$ and
$\min\left|\left|c^{G}\right|\right|$ denote the maximum and minimum nodal
concentration values from the ground truth result, respectively.

### Code and data availability.

The source code for our model and all input data are available for download
from a public software repository located at
https://github.com/truthlive/NeuronTransportLearning. All data generated
during this study can be reconstructed by running the source code.

## References

  * [1] Segev, I. & London, M.  Untangling dendrites with quantitative models.  _Science_ 290, 744–750 (2000).
  * [2] Swanger, S. A. & Bassell, G. J.  Dendritic protein synthesis in the normal and diseased brain.  _Neuroscience_ 232, 106–127 (2013).
  * [3] De Vos, K. J., Grierson, A. J., Ackerley, S. & Miller, C. C.  Role of axonal transport in neurodegenerative diseases.  _Nature Review Neuroscience_ 31, 151–173 (2008).
  * [4] Gunawardena, S. & Goldstein, L. S.  Polyglutamine diseases and transport problems: deadly traffic jams on neuronal highways.  _Archives of Neurology_ 62, 46–51 (2005).
  * [5] Millecamps, S. & Julien, J.-P.  Axonal transport deficits and neurodegenerative diseases.  _Nature Reviews Neuroscience_ 14, 161–176 (2013).
  * [6] Kononenko, N. L. _et al._ Retrograde transport of TrkB-containing autophagosomes via the adaptor AP-2 mediates neuronal complexity and prevents neurodegeneration.  _Nature Communications_ 8, 1–16 (2017).
  * [7] Zhang, H. _et al._ Modulation of AMPA receptor surface diffusion restores hippocampal plasticity and memory in Huntington’s disease models.  _Nature Communications_ 9, 1–16 (2018).
  * [8] Vale, R. D.  The molecular motor toolbox for intracellular transport.  _Cell_ 112, 467–480 (2003).
  * [9] Hirokawa, N., Niwa, S. & Tanaka, Y.  Molecular motors in neurons: transport mechanisms and roles in brain function, development, and disease.  _Neuron_ 68, 610–638 (2010).
  * [10] Franker, M. A. & Hoogenraad, C. C.  Microtubule-based transport–basic mechanisms, traffic rules and role in neurological pathogenesis.  _Journal of Cell Science_ 126, 2319–2329 (2013).
  * [11] Smith, D. & Simmons, R.  Models of motor-assisted transport of intracellular particles.  _Biophysical Journal_ 80, 45–68 (2001).
  * [12] Friedman, A. & Craciun, G.  A model of intracellular transport of particles in an axon.  _Journal of Mathematical Biology_ 51, 217–246 (2005).
  * [13] Craciun, G., Brown, A. & Friedman, A.  A dynamical system model of neurofilament transport in axons.  _Journal of Theoretical Biology_ 237, 316–322 (2005).
  * [14] Kuznetsov, A. & Avramenko, A.  A macroscopic model of traffic jams in axons.  _Mathematical Biosciences_ 218, 142–152 (2009).
  * [15] Kuznetsov, A., Avramenko, A. & Blinov, D.  Modeling the effect of a microtubule swirl on fast axonal transport.  _International Communications in Heat and Mass Transfer_ 37, 234–238 (2010).
  * [16] Hines, M. L. & Carnevale, N. T.  The NEURON simulation environment.  _Neural Computation_ 9, 1179–1209 (1997).
  * [17] Loew, L. M. & Schaff, J. C.  The Virtual Cell: a software environment for computational cell biology.  _TRENDS in Biotechnology_ 19, 401–406 (2001).
  * [18] Hughes, T., Cottrell, J. & Bazilevs, Y.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  _Computer Methods in Applied Mechanics and Engineering_ 194, 4135 – 4195 (2005).
  * [19] Li, A., Chai, X., Yang, G. & Zhang, Y. J.  An isogeometric analysis computational platform for material transport simulation in complex neurite networks.  _Molecular & Cellular Biomechanics_ 16, 123–140 (2019).
  * [20] Krizhevsky, A., Sutskever, I. & Hinton, G. E.  ImageNet classification with deep convolutional neural networks.  In _Advances in Neural Information Processing Systems_ , 1097–1105 (2012).
  * [21] LeCun, Y., Bengio, Y. & Hinton, G.  Deep learning.  _Nature_ 521, 436–444 (2015).
  * [22] Witten, I. H., Frank, E., Hall, M. A. & Pal, C. J.  _Data Mining: Practical Machine Learning Tools and Techniques_ (Morgan Kaufmann, 2016).
  * [23] Han, J., Jentzen, A. & Weinan, E.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ 115, 8505–8510 (2018).
  * [24] Raissi, M., Perdikaris, P. & Karniadakis, G. E.  Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ 378, 686–707 (2019).
  * [25] He, S. _et al._ Learning to predict the cosmological structure formation.  _Proceedings of the National Academy of Sciences_ 116, 13825–13832 (2019).
  * [26] Farimani, A. B., Gomes, J. & Pande, V. S.  Deep learning the physics of transport phenomena.  _arXiv Preprint:1709.02432_ (2017).
  * [27] Wiewel, S., Becher, M. & Thuerey, N.  Latent space physics: towards learning the temporal evolution of fluid flow.  _Computer Graphics Forum_ 38, 71–82 (2019).
  * [28] Li, A., Chen, R., Farimani, A. B. & Zhang, Y. J.  Reaction diffusion system prediction based on convolutional neural network.  _Scientific Reports_ 10, 1–9 (2020).
  * [29] Kipf, T. N. & Welling, M.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ (2016).
  * [30] Hamilton, W., Ying, Z. & Leskovec, J.  Inductive representation learning on large graphs.  In _Advances in Neural Information Processing Systems_ , 1024–1034 (2017).
  * [31] Monti, F. _et al._ Geometric deep learning on graphs and manifolds using mixture model CNNs.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 5115–5124 (2017).
  * [32] Fey, M., Lenssen, J. E., Weichert, F. & Müller, H.  SplineCNN: fast geometric deep learning with continuous B-spline kernels.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 869–877 (2018).
  * [33] Ogoke, F., Meidani, K., Hashemi, A. & Farimani, A. B.  Graph convolutional neural networks for body force prediction.  _arXiv preprint arXiv:2012.02232_ (2020).
  * [34] Karamad, M. _et al._ Orbital graph convolutional neural network for material property prediction.  _Physical Review Materials_ 4, 093801 (2020).
  * [35] Battaglia, P. W. _et al._ Relational inductive biases, deep learning, and graph networks.  _arXiv preprint arXiv:1806.01261_ (2018).
  * [36] Alet, F. _et al._ Graph element networks: adaptive, structured computation and memory.  _arXiv preprint arXiv:1904.09019_ (2019).
  * [37] Sanchez-Gonzalez, A. _et al._ Learning to simulate complex physics with graph networks.  _arXiv preprint arXiv:2002.09405_ (2020).
  * [38] Paszke, A. _et al._ Pytorch: an imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , 8024–8035 (Curran Associates, Inc., 2019).
  * [39] Fey, M. & Lenssen, J. E.  Fast graph representation learning with PyTorch Geometric.  _arXiv preprint arXiv:1903.02428_ (2019).
  * [40] Ascoli, G. A., Donohue, D. E. & Halavi, M.  Neuromorpho.org: a central resource for neuronal morphologies.  _Journal of Neuroscience_ 27, 9247–9251 (2007).
  * [41] Kingma, D. P. & Ba, J.  Adam: a method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * [42] Towns, J. _et al._ XSEDE: Accelerating scientific discovery.  _Computing in Science & Engineering_ 16, 62–74 (2014).
  * [43] Pan, S. J. & Yang, Q.  A survey on transfer learning.  _IEEE Transactions on Knowledge and Data Engineering_ 22, 1345–1359 (2009).
  * [44] Piegl, L. & Tiller, W.  _The NURBS Book_ (Springer Science & Business Media, 2012).
  * [45] Dimitri, R. _et al._ Isogeometric large deformation frictionless contact using T-splines.  _Computer Methods in Applied Mechanics and Engineering_ 269, 394–414 (2014).
  * [46] Benson, D., Bazilevs, Y., Hsu, M.-C. & Hughes, T. J.  Isogeometric shell analysis: the Reissner–Mindlin shell.  _Computer Methods in Applied Mechanics and Engineering_ 199, 276–289 (2010).
  * [47] Casquero, H. _et al._ Seamless integration of design and Kirchhoff–Love shell analysis using analysis-suitable unstructured T-splines.  _Computer Methods in Applied Mechanics and Engineering_ 360, 112765 (2020).
  * [48] Casquero, H. _et al._ Arbitrary-degree T-splines for isogeometric analysis of fully nonlinear Kirchhoff–Love shells.  _Computer-Aided Design_ 82, 140–153 (2017).
  * [49] Zhang, Y., Bazilevs, Y., Goswami, S., Bajaj, C. L. & Hughes, T. J.  Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow.  _Computer Methods in Applied Mechanics and Engineering_ 196, 2943–2959 (2007).
  * [50] Zhang, Y. _et al._ An atlas-based geometry pipeline for cardiac Hermite model construction and diffusion tensor reorientation.  _Medical Image Analysis_ 16, 1130–1141 (2012).
  * [51] Urick, B., Sanders, T. M., Hossain, S. S., Zhang, Y. J. & Hughes, T. J.  Review of patient-specific vascular modeling: template-based isogeometric framework and the case for CAD.  _Archives of Computational Methods in Engineering_ 26, 381–404 (2019).
  * [52] Yu, Y., Zhang, Y. J., Takizawa, K., Tezduyar, T. E. & Sasaki, T.  Anatomically realistic lumen motion representation in patient-specific space–time isogeometric flow analysis of coronary arteries with time-dependent medical-image data.  _Computational Mechanics_ 65, 395–404 (2020).
  * [53] Zhang, Y. J.  _Geometric Modeling and Mesh Generation from Scanned Images_ , vol. 6 (CRC Press, 2016).
  * [54] Bazilevs, Y., Calo, V. M., Zhang, Y. & Hughes, T. J.  Isogeometric fluid–structure interaction analysis with applications to arterial blood flow.  _Computational Mechanics_ 38, 310–322 (2006).
  * [55] Casquero, H., Liu, L., Bona-Casas, C., Zhang, Y. & Gomez, H.  A hybrid variational-collocation immersed method for fluid-structure interaction using unstructured T-splines.  _International Journal for Numerical Methods in Engineering_ 105, 855–880 (2016).
  * [56] Casquero, H., Zhang, Y. J., Bona-Casas, C., Dalcin, L. & Gomez, H.  Non-body-fitted fluid–structure interaction: divergence-conforming B-splines, fully-implicit dynamics, and variational formulation.  _Journal of Computational Physics_ 374, 625–653 (2018).
  * [57] Anitescu, C., Jia, Y., Zhang, Y. J. & Rabczuk, T.  An isogeometric collocation method using superconvergent points.  _Computer Methods in Applied Mechanics and Engineering_ 284, 1073–1097 (2015).
  * [58] Casquero, H., Liu, L., Zhang, Y., Reali, A. & Gomez, H.  Isogeometric collocation using analysis-suitable T-splines of arbitrary degree.  _Computer Methods in Applied Mechanics and Engineering_ 301, 164–186 (2016).
  * [59] Wei, X., Zhang, Y., Liu, L. & Hughes, T. J.  Truncated T-splines: fundamentals and methods.  _Computer Methods in Applied Mechanics and Engineering_ 316, 349–372 (2017).
  * [60] Wei, X., Zhang, Y. J. & Hughes, T. J.  Truncated hierarchical tricubic $C^{0}$ spline construction on unstructured hexahedral meshes for isogeometric analysis applications.  _Computers & Mathematics with Applications_ 74, 2203–2220 (2017).
  * [61] Wei, X. _et al._ Blended B-spline construction on unstructured quadrilateral and hexahedral meshes with optimal convergence rates in isogeometric analysis.  _Computer Methods in Applied Mechanics and Engineering_ 341, 609–639 (2018).
  * [62] Li, X., Wei, X. & Zhang, Y. J.  Hybrid non-uniform recursive subdivision with improved convergence rates.  _Computer Methods in Applied Mechanics and Engineering_ 352, 606–624 (2019).
  * [63] Scarselli, F., Gori, M., Tsoi, A. C., Hagenbuchner, M. & Monfardini, G.  The graph neural network model.  _IEEE Transactions on Neural Networks_ 20, 61–80 (2008).
  * [64] Zhang, M. & Chen, Y.  Link prediction based on graph neural networks.  In _Advances in Neural Information Processing Systems_ , 5165–5175 (2018).
  * [65] Geng, X. _et al._ Spatiotemporal multi-graph convolution network for ride-hailing demand forecasting.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, 3656–3663 (2019).
  * [66] Battaglia, P., Pascanu, R., Lai, M., Rezende, D. J. _et al._ Interaction networks for learning about objects, relations and physics.  In _Advances in Neural Information Processing Systems_ , 4502–4510 (2016).
  * [67] Santoro, A. _et al._ A simple neural network module for relational reasoning.  In _Advances in Neural Information Processing Systems_ , 4967–4976 (2017).
  * [68] Khamsi, M. A. & Kirk, W. A.  _An Introduction to Metric Spaces and Fixed Point Theory_ , vol. 53 (John Wiley & Sons, 2011).
  * [69] Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O. & Dahl, G. E.  Neural message passing for quantum chemistry.  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , 1263–1272 (2017).
  * [70] Wang, X., Girshick, R., Gupta, A. & He, K.  Non-local neural networks.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 7794–7803 (2018).

## Acknowledgements

The authors acknowledge the support of NSF grants CMMI-1953323 and
CBET-1804929. This work used the Extreme Science and Engineering Discovery
Environment (XSEDE), which is supported by the NSF grant ACI-1548562.
Specifically, it used the Bridges system, which is supported by the NSF grant
ACI-1445606, at the Pittsburgh Supercomputing Center (PSC).

## Author contributions statement

A.L. designed the model and carried out the calculations and analysis with the
instruction from Y.J.Z.. A.B.F. provided advice on the selection of machine
learning methods. All authors contributed to the writing, discussions and
revisions of the manuscript.

## Additional information

### Competing interests.

The authors declare no competing interests.

## Supplementary Figures

Figure S1: The graph extraction of the pipe and bifurcation structures. (A, B)
The input control meshes of the pipe and the bifurcation. (C) The extraction
template for circular cross sections with 17 extracted nodes labeled in red.
(D) The extraction template for branch cross sections with 23 extracted nodes
labeled in red. (E, F) The extracted graphs of the pipe and the bifurcation.
Figure S2: The test loss vs epoch curves for two GNN simulators. 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 all-hexahedral
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.