# Electronic correlations and Fermi liquid behavior of intermediate-band
states in titanium-doped silicon

A. Östlin  Theoretical Physics III, Center for Electronic Correlations and
Magnetism, Institute of Physics, University of Augsburg, D-86135 Augsburg,
Germany  Augsburg Center for Innovative Technologies, University of Augsburg,
D-86135 Augsburg, Germany  L. Chioncel  Augsburg Center for Innovative
Technologies, University of Augsburg, D-86135 Augsburg, Germany  Theoretical
Physics III, Center for Electronic Correlations and Magnetism, Institute of
Physics, University of Augsburg, D-86135 Augsburg, Germany

###### Abstract

We study the nature of the electronic states in the intermediate band formed
by interstitial titanium in silicon. Our single-site description combines
effects of electronic correlations, captured by dynamical mean-field theory,
and disorder, modeled using the coherent potential approximation and the
typical medium mean-field theory. For all studied concentrations an extended
metallic state with a strongly depleted density of states at the Fermi level
is obtained. The self-energy is characteristic to Fermi-liquids and for
certain temperatures reveals the existence of coherent quasi-particles.

The metal-to-insulator transition (MIT) is a rich field of research within
condensed-matter physics, where many different underlying mechanisms (disorder
and/or interaction) has been proposed [1, 2]. Anderson [3] introduced a
mechanism which now bears his name, while studying doping of phosphor in
silicon, where the transition is due to the randomness of the impurities. At
about the same time, Mott [4] proposed that screening effects would become
more important as the number of impurities increases, leading to a MIT. A
third possibility is the Mott-Hubbard transition [1], in which an on-site
intra-atomic electron-electron Coulomb interaction competes with inter-atomic
hopping, leading potentially to an insulating state. In principle, all three
of these mechanisms are present in doped semiconductors, and to disentangle
one from the other is a non-trivial task.

From the family of doped semiconductors, one promising candidate for the next
generation of solar cells are the intermediate-band (IB) photovoltaics [5, 6,
7]. In this class of materials, deep-level impurities are used to create a
narrow band situated inside the gap of a semiconductor host. The IB allows
two-photon absorption processes, which are beneficial for the conversion
efficiency of solar-cell devices. However, the deep-level impurity may also
give rise to a localized charge around itself, which in turn distorts the
lattice. This can lead to non-radiative recombination processes, as discussed
by Shockley, Reed [8], and Hall [9], where the efficiency is reduced due to
phonons. It has been argued that the non-radiative recombination can be
counteracted by increasing the doping of impurities to the point where the
states in the IB becomes delocalized, making the charge spread out [10, 11].
Therefore, it is an important question to answer if the IB states are
localized or delocalized, or in other words, if the system is insulating or
metallic.

An often studied IB photovoltaic is silicon doped with titanium [12, 13, 14].
During the last two decades, there has been a substantial research interest in
this system, with conflicting results and interpretations. It seems that both
the crystal growth techniques used to grow the silicon host, and the
implantation method for the Ti, play important roles. Both experiment and
theory indicates that Ti impurities predominantly occupy interstitial
positions in the diamond structure of Si [11, 15, 16, 17]. Since the
solubility limit of Ti in Si is relatively low [18, 19], out-of-equilibrium
implantation techniques have to be used to achieve so-called hyperdoping. The
Mott limit, where the states in the IB would turn from localized to
delocalized, is estimated to be $N_{c}=6\times 10^{19}$ cm-3 [11], therefore
hyperdoping is necessary to reach the MIT. High Ti concentrations have been
reported, some showing as high as $10^{19}$ and $10^{21}$ cm-3 [12, 20].
However, even though good conductivity and sub-bandgap photoresponse has been
demonstrated, the specific detectivity was found to be rather low [14]. Liu
_et al._ [16] used two different Ti-implantation techniques, where one sample
was completely insulating for all concentrations, while the other had a
decreasing sheet resistance with increasing Ti concentration. In the latter
case, the decrease in conductivity was attributed to percolation in Ti-rich
cellular walls, and not to the MIT. This brings into question the degree of
spatial homogeneity of Ti one can expect after hyperdoping. Recent studies,
using ion implantation and pulsed-laser melting, came to the conclusion that
the possible concentrations are well below the Mott limit [17, 21].

On the theoretical side, in an attempt to formulate a mean-field theory of
Anderson localization involving one-particle quantities, it was realized that
the typical density of states (TDOS) vanishes as the disorder strength
increases, and therefore can serve as a possible order parameter for
localization [22, 23, 24]. Such an effect can not be captured by using the
arithmetic average of the self-consistent coherent potential approximation
(CPA) [25, 26]. Within the proposed typical medium theory (TMT) [23],
geometric averaging over disorder was used, instead of the arithmetic average.
A combined density functional theory (DFT) [27] calculation followed by the
cluster extension of the TMT, the so-called typical medium dynamical cluster
approximation TM-DCA [28, 29] recently reported that in the concentration
range of $1.0\times 10^{20}$ to $5.0\times 10^{19}$ cm-3 Ti impurities in Si
[30] the IB electronic states within the gap undergo a localized to extended
transition.

As the incomplete $3d$ shell of Ti is subject to important electronic
correlation effects, we investigate here the competition of electron-electron
interaction and localization of IB electronic states in Ti-doped Si. On the
methodological side we report the combined implementation of TMT and dynamical
mean-field theory (DMFT) [31, 32], to describe electronic correlation, into
the DFT framework. This is a natural extension of a previous single-site TMT
formulation using the exact muffin-tin orbitals (EMTO) basis set [33], which
however departs from the original single-site TMT formulation as will be
explained below. We compare spectral functions of DFT(CPA/TMT) [33], with and
without electronic correlation [34], with our present TMT+DMFT implementation.
According to the present results the system remains metallic for the
experimentally relevant Hubbard $U$ and $J$ parameters with a significantly
reduced density of states at the Fermi level and a non-zero but reduced TMT
order parameter.

Figure 1: CPA/TMT+DMFT flow diagram. The path operator $\tilde{g}$ for the
effective medium is arithmetically averaged in single-site CPA, while it is
geometrically averaged in single-site TMT. The potential parameters $P$ and
structure constants $S$ parametrize the Green’s function [34]. The full charge
self-consistency is achieved in the DFT external loop: the charge density
$n({\bf r})$ is passed into the combined disorder and many-body solver, while
the self-energy $\Sigma$ is stored for the next many-body iteration, and the
correction $\Delta n({\bf r})$ to the real space charge density is returned
into the outside DFT loop.

In order to clarify the implementation let us state the key idea of the TMT
namely that the typical values of the random quantities should be associated
with physical observables [3]. The typical value of a variable $X$ is
represented by its geometric mean $\langle
X\rangle_{\mathrm{g}}\equiv\mathrm{exp}\
[\langle\mathrm{ln}X(\varepsilon)\rangle_{\mathrm{a}}]$, where
$X(\varepsilon)$ is a function of the random on-site energy $\varepsilon$ and
by $\langle...\rangle_{\mathrm{a}}$ we denote the arithmetic average. The
typical density of states reads: $\rho_{g}(E)$=$\mathrm{exp}\
[\langle\ln\rho(E,\varepsilon)\rangle_{\mathrm{a}}]$. The typical value of the
local Green’s function is obtained as a Hilbert transform:
$G_{g}(z)=\int_{-\infty}^{\infty}\rho_{g}(E)/(z-E)\mathrm{d}E$, which in many
situations lead to numerical difficulties. In our TMT formulation we
geometrically average the full scattering path operator ($\tilde{g}$, see Fig.
1) in such a way we avoid the numerical instabilities associated with the
Hilbert transform. We follow the implementation of Ref. [34], which
parametrize the Green’s function (in terms of potential parameters $P$ and
structure constants $S$) and thus does not require any analytical continuation
within the DFT+DMFT self-consistency loops [35]. The schematic flow diagram is
shown in Fig. 1 and follows that presented in Ref. [34].

In order to solve the DMFT impurity problem, the perturbative spin-polarized
$T$-matrix fluctuation-exchange (SPTFLEX) solver [36, 37, 38, 39, 40] is used.
In this solver, the electron-electron interaction term can be considered in a
full spin and orbital rotationally invariant form, viz.
$\frac{1}{2}\sum_{{i\\{m,\sigma\\}}}U_{mm^{\prime}m^{\prime\prime}m^{\prime\prime\prime}}c^{{\dagger}}_{im\sigma}c^{{\dagger}}_{im^{\prime}\sigma^{\prime}}c_{im^{\prime\prime\prime}\sigma^{\prime}}c_{im^{\prime\prime}\sigma}$.
Here, $c_{im\sigma}/c^{\dagger}_{im\sigma}$ annihilates/creates an electron
with spin $\sigma$ on the orbital $m$ at the lattice site $i$. The Coulomb
matrix elements $U_{mm^{\prime}m^{\prime\prime}m^{\prime\prime\prime}}$ are
expressed in the usual way in terms of Slater integrals [1]. The
simplifications of the computational procedure in reformulating the FLEX [36]
as a DMFT impurity solver [37, 38, 39, 41] consists in neglecting dynamical
interactions in the particle-particle channel, considering only static (of
$T-$matrix type) renormalization of the effective interactions. As some
interaction effects are already included in the DFT exchange-correlation
functional, the so-called “double-counted” terms must be subtracted. To
achieve this, we employ the fully-localized limit double-counting scheme [42].
The Matsubara frequencies $i\omega_{n}=(2n+1)i\pi T$, with $n=0,1,2,\cdots$,
were truncated after 1024 frequencies, and different values for the
temperature $T$ were considered. The densities of state were computed along a
horizontal contour shifted away from the real-energy axis by about $10$ meV.
At the end of the self-consistent DFT calculations, the self-energy
$\Sigma(i\omega_{n})$ was analytically continued by a Padé approximant [43,
44] to the horizontal contour.

Silicon crystallizes in the diamond structure (space group 227,
Fd$\overline{3}$m), with a lattice constant $a=5.4309$ Å. The Si atoms occupy
the $(0,0,0)$ and $a(\frac{1}{4},\frac{1}{4},\frac{1}{4})$ positions, while
the $a(\frac{1}{2},\frac{1}{2},\frac{1}{2})$ and
$a(\frac{3}{4},\frac{3}{4},\frac{3}{4})$ are unoccupied. We denote the
unoccupied positions as Em (for “empty”) and consider the Ti-interstitial
impurities as randomly placed at the $a(\frac{1}{2},\frac{1}{2},\frac{1}{2})$
lattice-sites. We view the Ti and the Em as forming a “binary alloy”,
occupying the $a(\frac{1}{2},\frac{1}{2},\frac{1}{2})$-site with the formula
TicEm1-c, where $0\leq c\leq 1$ denotes the concentration of Ti. Throughout
the paper, unless stated otherwise, we use $c=0.001$, since this is a
representative concentration close to the predicted Mott limit [11]. We use
the experimental lattice constant throughout all calculations, since for the
investigated concentrations of Ti impurities the change in unit cell volume
upon doping is negligible [15, 45]. The single-site method that we use here
can not, in general, capture the effect of local lattice relaxations. However,
in previous studies the position of the IB was found to be fairly unaffected
by the local lattice relaxation due to the Ti-impurities [15, 30]. The Green’s
function is computed for 16 energy points along a semicircular contour with a
diameter of 1 Ry, with the basis consisting of $s$-, $p$-, and $d$-orbitals.
The DFT convergence is achieved already at about $5\cdot 10^{3}$ k-points
however DOS plots were produced using an integration over the irreducible
Brillouin zone a mesh consisting of more than $10^{4}$ k-points. The
generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE)
parametrization [46] of the exchange-correlation functional was used
throughout this paper. The GGA was used by Sánchez _et al._ [15] in studying
Ti impurities in Si using supercells.

Figure 2: Concentration weighted Si partial DOS using CPA (black line) and TMT
(orange line). Inset: Concentration-weighted partial DOS for interstitial Ti,
at a concentration of $0.1\%$. The partial density of states of pure Si is
shown as a shaded area.

In Fig. 2 we show the partial LDOS for Si and Ti (inset) in the non-
interacting case using the CPA and TMT configurational averages. The IB states
are formed by the hybridization of Ti-$3d$ states and the Si-$p$ states. By
orbital decomposition we confirm that at the Fermi level $E_{F}$, Ti-$t_{2g}$
orbitals are predominant. The impurity band in the CPA overlaps with the
bottom conduction band, which is a consequence of the poor performance of the
traditional exchange correlation functional, which leads to an underestimated
gap. In general the arithmetic average is larger than the geometric one [47],
however for such a small concentration the average LDOS for the CPA and TMT
are quite similar. This might be attributed to the fact that Ti impurities are
subject to similar hybridization with the tails of the Si-atoms on the vacated
interstitial site. We investigate further the behaviour of the TDOS, the so-
called order parameter which monitors the presence of precursors for Anderson
localization [48, 23, 24]. In the absence of electronic correlations this
quantity is computed according to:
$\rho_{g}(E)=\rho^{c}_{\mathrm{Ti}}(E)\rho^{1-c}_{\mathrm{Em}}(E)$, and
represents the geometrical average of the Ti-LDOS and the LDOS produced by all
other atoms (mainly Si) at the position of the empty-site. Its magnitude
remains essentially unmodified $\rho_{g}(E_{F})\approx 0.008$ (eV)-1 in the
entire range of studied concentrations, and thus its non-zero value signals
the lack of a disorder-driven MIT.

The effects of the mean-field decoupling of the Hubbard interaction has been
addressed using the DFT+$U$ method [49]. A Hubbard correction to the $3d$
orbitals in the range of $U\approx 3.0-5.0$ eV give a good agreement for the
magnitude of the migration barriers and the position of donor levels for
certain transition metal impurities [49, 45] including Ti. Given the spread in
the magnitude of the $U$ parameter we have performed (CPA/TMT)+DMFT
computations for several $U$ and $J$.

Figure 3: Ti partial density of states, for a doping level of $c=0.1\%$.
CPA+DMFT (black lines) and TMT+DMFT (orange lines) for different values of the
Coulomb $U$ and $J$ parameters for $T=400$ K.

In Fig. 3 we show the Ti partial LDOS for different averaging schemes, and for
different $U$ and $J$ parameters, at $T=400$ K. Electronic correlation effects
on the Ti-$3d$ states lead to a spectral weight transfer towards the bottom of
the conduction band (towards higher energies). Most significantly, depending
on the strength of Hubbard parameters a major reduction of about two orders of
magnitude of the Ti-LDOS at the Fermi level is obtained in the range of
$\rho_{\mathrm{Ti}}(E_{F})\approx 10^{-3}$ to $10^{-5}$ states/eV.
Nevertheless the total DOS, $\rho(E_{F})$, is not zero but with an absolute
magnitude of about $10^{-2}$ states/eV. The CPA+DMFT methodology allows also
to study the temperature dependence of the partial LDOS at the impurity site.
For $U=5.0$ eV, $J=1.0$ eV we found that in the vicinity of the Fermi level
Ti-LDOS decreases with temperature, however its magnitude is significantly
smaller than the corresponding Em-LDOS (not shown). Overall no significant
temperature dependence of the total DOS at $E_{F}$ is obtained.

Figure 4: Temperature dependence of the orbital resolved Im
$\Sigma_{e_{g}/t_{2g}}(i\omega_{n})$ within the CPA+DMFT and at $U=5.0$ eV,
J=1.0 eV. Squares/circles denote the $e_{g}/t_{2g}$-orbitals character. Lines
are guides to the eye. Inset: Im $\Sigma_{e_{g}/t_{2g}}(i\omega_{0})$ plotted
versus temperature. The dashed line is a linear fit to the data.

In metallic alloys the self-energy carries information about the relaxation
processes. In particular the imaginary part of the DMFT self-energy encodes
the Coulomb scattering effects. In disordered systems in the absence of
interactions the self-energy takes into account the fluctuations in the on-
site energies beyond the effective medium fields of CPA and TMT, and carries
the information about the scattering on impurities. In the case of local Fermi
liquids it was recently shown that even in the presence of a constant
scattering rate (which may mimic impurity scattering) into the DMFT self-
energy the system still displays a Fermi-liquid behavior [50].

In Fig. 4 we show the Ti $3d$-orbital-resolved self-energies as functions of
the Matsubara frequencies from the CPA+DMFT computation. Similar results were
obtained from TMT+DMFT. The Fermi liquid state is characterized by a linear
dependence of the imaginary part of the self-energy on imaginary frequencies
$i\omega_{n}$: $\mathrm{Im}\Sigma(i\omega_{n})=-\Gamma-(Z^{-1}-1)\omega_{n}$,
where $Z$ is the quasiparticle spectral weight and $\Gamma$ the quasiparticle
damping (inverse quasiparticle life time). The
$\mathrm{Im}\Sigma_{e_{g}}(i\omega_{n})$ approaches zero linearly for all
temperatures in the region of low Matsubara frequencies, thus the
$e_{g}$-electrons are long-lived quasiparticles, i.e. their scattering rate
vanishes as the Fermi surface is approached. A downturn in the
$\mathrm{Im}\Sigma_{t_{2g}}(i\omega_{n})$ is seen in the temperature range
from 1000 K down to 300 K. However, for temperatures between
$200\,\mathrm{K}\leq T\leq 300$ K, the
$\mathrm{Im}\Sigma_{t_{2g}}(i\omega_{n})$ is changing its curvature, so one
expects that in this temperature range a cross over to the Fermi liquid regime
occurs. The inset of Fig. 4 shows that Im $\Sigma_{e_{g}/t_{2g}}(i\omega_{0})$
versus temperature can be fitted to a line. In view of the “first Matsubara
frequency” rule [51] such a linear dependence provides an estimate for the
crossover temperature below which quasiparticles are coherent. While
$e_{g}$-electrons are coherent quasiparticles in the entire temperature range
the $t_{2g}$-electrons behave as coherent quasiparticles only below $\approx
200$ K. These results provide a clear indication of the Fermi-liquid character
of the electronic system. The quasiparticle mass enhancement $m^{\star}/m$,
with respect to the band mass $m$, can be approximately computed from the
Matsubara frequencies in the zero-temperature limit as
$m^{\star}/m=1-\frac{\mathrm{Im}\Sigma(i\omega_{n})}{\omega_{n}}|_{\omega_{n}\rightarrow
0}$. For the Hubbard parameter $U=5$ eV, we estimate a mass enhancement of
$m^{\star}/m\approx 8$ for the $t_{2g}$ orbitals at $150$ K, which confirms
the presence of important correlation effects.

In summary, we modeled randomness and local Coulomb interaction effects in
titanium-doped silicon. We compare results of CPA+DMFT [34] with the present
TMT+DMFT implementation and address the controversy whether the induced IB
states, formed as consequence of Ti-$d$($t_{2g}$) and Si-$p$ orbital
hybridization, are localized or extended. We demonstrate that within the
single-site version of TMT-EMTO, in the absence of electronic correlations,
the order parameter that describes the precursors of Anderson localization is
non-zero, so that no transition can be captured. In the presence of
interactions we report a significant reduction of the one-particle spectra at
$E_{F}$, however the electronic states remain coherent quasiparticles at least
up to $200$ K, despite a significantly large Hubbard $U$ parameter. The
drastic reduction of $\rho(E_{F})$ is a predominant correlation effect being
qualitatively and even quantitatively similar between CPA+DMFT and TMT+DMFT.
Thus, if a metal to insulator transition exists this should be triggered by
effects which are beyond the combined single-site disorder and correlation
considered here. Note that recent experimental findings, points towards
granular and percolation physics being of importance in this system, which
requires modeling beyond local theories.

We have benefited from discussions with H. Terletska, Y. Zhang and L. Vitos.
Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) — Grant number 107745057 — TRR 80.

## References

  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998).
  * Dobrosavljevic  _et al._ [2012] V. Dobrosavljevic, N. Trivedi, and J. M. Valles Jr., eds.,  _Conductor-Insulator Quantum Phase Transitions_ (Oxford University Press, Oxford, 2012).
  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).
  * Mott [1968] N. F. Mott, Metal-insulator transition, Rev. Mod. Phys. 40, 677 (1968).
  * Luque and Martí [1997] A. Luque and A. Martí, Increasing the efficiency of ideal solar cells by photon induced transitions at intermediate levels, Phys. Rev. Lett. 78, 5014 (1997).
  * Luque and Martí [2010] A. Luque and A. Martí, The intermediate band solar cell: Progress toward the realization of an attractive concept, Advanced Materials 22, 160 (2010).
  * Luque  _et al._ [2012] A. Luque, A. Martí, and C. Stanley, Understanding intermediate-band solar cells, Nature Photonics 6, 146 (2012).
  * Shockley and Read [1952] W. Shockley and W. T. Read, Statistics of the recombinations of holes and electrons, Phys. Rev. 87, 835 (1952).
  * Hall [1952] R. N. Hall, Electron-hole recombination in germanium, Phys. Rev. 87, 387 (1952).
  * Lang and Henry [1975] D. V. Lang and C. H. Henry, Nonradiative Recombination at Deep Levels in GaAs and GaP by Lattice-Relaxation Multiphonon Emission, Phys. Rev. Lett. 35, 1525 (1975).
  * Luque  _et al._ [2006] A. Luque, A. Martí, E. Antolín, and C. Tablero, Intermediate bands versus levels in non-radiative recombination, Physica B: Condensed Matter 382, 320 (2006).
  * Antolín  _et al._ [2009] E. Antolín, A. Martí, J. Olea, D. Pastor, G. González-Díaz, I. Mártil, and A. Luque, Lifetime recovery in ultrahighly titanium-doped silicon for the implementation of an intermediate band material, Applied Physics Letters 94, 042115 (2009), https://doi.org/10.1063/1.3077202 .
  * Olea  _et al._ [2016] J. Olea, E. López, E. Antolín, A. Martí, A. Luque, E. García-Hemme, D. Pastor, R. García-Hernansanz, A. del Prado, and G. González-Díaz, Room temperature photo-response of titanium supersaturated silicon at energies over the bandgap, J. Phys. D: Appl. Phys. 49, 055103 (2016).
  * Wang and Berencén [2021] M. Wang and Y. Berencén, Room-temperature infrared photoresponse from ion beam–hyperdoped silicon, physica status solidi (a) 218, 2000260 (2021).
  * Sánchez  _et al._ [2009] K. Sánchez, I. Aguilera, P. Palacios, and P. Wahnón, Assessment through first-principles calculations of an intermediate-band photovoltaic material based on Ti-implanted silicon: Interstitial versus substitutional origin, Phys. Rev. B 79, 165203 (2009).
  * Liu  _et al._ [2018] F. Liu, M. Wang, Y. Berencén, S. Prucnal, M. Engler, R. Hübner, Y. Yuan, R. Heller, R. Böttger, L. Rebohle, W. Skorupa, M. Helm, and S. Zhou, On the insulator-to-metal transition in titanium-implanted silicon, Scientific Reports 8, 4164 (2018).
  * Lim  _et al._ [2021] S. Q. Lim, A. J. Akey, E. Napolitani, P. K. Chow, J. M. Warrender, and J. S. Williams, A critical evaluation of Ag- and Ti-hyperdoped Si for Si-based infrared light detection, Journal of Applied Physics 129, 065701 (2021), https://doi.org/10.1063/5.0035620 .
  * Hocine and Mathiot [1988] S. Hocine and D. Mathiot, Titanium diffusion in silicon, Applied Physics Letters 53, 1269 (1988), https://doi.org/10.1063/1.100446 .
  * Mathiot and Barbier [1991] D. Mathiot and D. Barbier, Solubility enhancement of metallic impurities in silicon by rapid thermal annealing, Journal of Applied Physics 69, 3878 (1991), https://doi.org/10.1063/1.348444 .
  * Olea  _et al._ [2012] J. Olea, D. Pastor, E. García-Hemme, R. García-Hernansanz, Álvaro del Prado, I. Mártil, and G. González-Díaz, Low temperature intermediate band metallic behavior in Ti implanted Si, Thin Solid Films 520, 6614 (2012).
  * Akey  _et al._ [2021] A. J. Akey, J. Mathews, and J. M. Warrender, Maximum Ti concentrations in Si quantified with atom probe tomography (APT), Journal of Applied Physics 129, 175701 (2021), https://doi.org/10.1063/5.0029981 .
  * Dobrosavljević and Kotliar [1997] V. Dobrosavljević and G. Kotliar, Mean field theory of the Mott-Anderson transition, Phys. Rev. Lett. 78, 3943 (1997).
  * Dobrosavljević  _et al._ [2003] V. Dobrosavljević, A. A. Pastor, and B. K. Nikolić, Typical medium theory of Anderson localization: A local order parameter approach to strong-disorder effects, EPL 62, 76 (2003).
  * Dobrosavljević [2010] V. Dobrosavljević, Typical medium theory of Mott-Anderson localization, Int. J. Mod. Phys. B 24, 1680 (2010).
  * Taylor [1967] D. Taylor, Vibrational Properties of Imperfect Crystals with Large Defect Concentrations, Phys. Rev. 156, 1017 (1967).
  * Soven [1967] P. Soven, Coherent-Potential Model of Substitutional Disordered Alloys, Phys. Rev. 156, 809 (1967).
  * Jones [2015] R. O. Jones, Density functional theory: Its origins, rise to prominence, and future, Rev. Mod. Phys. 87, 897 (2015).
  * Ekuma  _et al._ [2014] C. E. Ekuma, H. Terletska, K.-M. Tam, Z.-Y. Meng, J. Moreno, and M. Jarrell, Typical medium dynamical cluster approximation for the study of Anderson localization in three dimensions, Phys. Rev. B 89, 081107 (2014).
  * Zhang  _et al._ [2015] Y. Zhang, H. Terletska, C. Moore, C. Ekuma, K.-M. Tam, T. Berlijn, W. Ku, J. Moreno, and M. Jarrell, Study of multiband disordered systems using the typical medium dynamical cluster approximation, Phys. Rev. B 92, 205111 (2015).
  * Zhang  _et al._ [2018] Y. Zhang, R. Nelson, K.-M. Tam, W. Ku, U. Yu, N. S. Vidhyadhiraja, H. Terletska, J. Moreno, M. Jarrell, and T. Berlijn, Origin of localization in Ti-doped Si, Phys. Rev. B 98, 174204 (2018).
  * Metzner and Vollhardt [1989] W. Metzner and D. Vollhardt, Correlated Lattice Fermions in d=$\infty$ dimensions, Phys. Rev. Lett. 62, 324 (1989).
  * Kotliar and Vollhardt [2004] G. Kotliar and D. Vollhardt, Strongly correlated materials: Insights from dynamical mean-field theory, Phys. Today 57, 53 (2004).
  * Östlin  _et al._ [2020] A. Östlin, Y. Zhang, H. Terletska, F. Beiuşeanu, V. Popescu, K. Byczuk, L. Vitos, M. Jarrell, D. Vollhardt, and L. Chioncel, Ab initio typical medium theory of substitutional disorder, Phys. Rev. B 101, 014210 (2020).
  * Östlin  _et al._ [2018] A. Östlin, L. Vitos, and L. Chioncel, Correlated electronic structure with uncorrelated disorder, Phys. Rev. B 98, 235135 (2018).
  * Östlin  _et al._ [2017] A. Östlin, L. Vitos, and L. Chioncel, Analytic continuation-free Green’s function approach to correlated electronic structure calculations, Phys. Rev. B 96, 125156 (2017).
  * Bickers and Scalapino [1989] N. E. Bickers and D. J. Scalapino, Conserving approximation for strongly fluctuating electron systems. I. formalism and calculational approach, Ann. Phys. (N. Y.) 193, 206 (1989).
  * Lichtenstein and Katsnelson [1998] A. I. Lichtenstein and M. I. Katsnelson, Ab initio calculations of quasiparticle band structure in correlated systems: LDA++ approach, Phys. Rev. B 57, 6884 (1998).
  * Katsnelson and Lichtenstein [1999] M. I. Katsnelson and A. I. Lichtenstein, LDA++ approach to the electronic structure of magnets: correlation effects in iron, J. Phys.: Condens. Matter 11, 1037 (1999).
  * Pourovskii  _et al._ [2005] L. V. Pourovskii, M. I. Katsnelson, and A. I. Lichtenstein, Correlation effects in electronic structure of actinide monochalcogenides, Phys. Rev. B 72, 115106 (2005).
  * Grånäs  _et al._ [2012] O. Grånäs, I. Di Marco, P. Thunström, L. Nordström, O. Eriksson, T. Björkman, and J. Wills, Charge self-consistent dynamical mean-field theory based on the full-potential linear muffin-tin orbital method: Methodology and applications, Computational Materials Science 55, 295–302 (2012).
  * Katsnelson  _et al._ [2008] M. I. Katsnelson, V. Y. Irkhin, L. Chioncel, A. I. Lichtenstein, and R. A. de Groot, Half-metallic ferromagnets: From band structure to many-body effects, Rev. Mod. Phys. 80, 315 (2008).
  * Petukhov  _et al._ [2003] A. G. Petukhov, I. I. Mazin, L. Chioncel, and A. I. Lichtenstein, Correlated metals and the $LDA+U$ method, Phys. Rev. B 67, 153106 (2003).
  * Vidberg and Serene [1977] H. J. Vidberg and J. W. Serene, Solving the Eliashberg equations by means of N-point Padé approximants, J. Low Temp. Phys. 29, 179 (1977).
  * Östlin  _et al._ [2012] A. Östlin, L. Chioncel, and L. Vitos, One-particle spectral function and analytic continuation for many-body implementation in the exact muffin-tin orbitals method, Phys. Rev. B 86, 235107 (2012).
  * Marinopoulos  _et al._ [2015] A. G. Marinopoulos, P. Santos, and J. Coutinho, DFT+${U}$ study of electrical levels and migration barriers of early $3d$ and $4d$ transition metals in silicon, Phys. Rev. B 92, 075124 (2015).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Bullen [2003] P. S. Bullen,  _Handbook of Means and Their Inequalities_ (Kluwer Academic Publishers, Dordrecht, 2003).
  * Schubert  _et al._ [2010] G. Schubert, J. Schleede, K. Byczuk, H. Fehske, and D. Vollhardt, Distribution of the local density of states as a criterion for Anderson localization: Numerically exact results for various lattices in two and three dimensions, Phys. Rev. B 81, 155106 (2010).
  * Markevich  _et al._ [2014] V. P. Markevich, S. Leonard, A. R. Peaker, B. Hamilton, A. G. Marinopoulos, and J. Coutinho, Titanium in silicon: Lattice positions and electronic properties, Appl. Phys. Lett. 104, 152105 (2014), https://doi.org/10.1063/1.4871702 .
  * Berthod  _et al._ [2013] C. Berthod, J. Mravlje, X. Deng, R. Žitko, D. van der Marel, and A. Georges, Non-Drude universal scaling laws for the optical response of local Fermi liquids, Phys. Rev. B 87, 115109 (2013).
  * Chubukov and Maslov [2012] A. V. Chubukov and D. L. Maslov, First-Matsubara-frequency rule in a Fermi liquid. I. Fermionic self-energy, Phys. Rev. B 86, 155136 (2012).


# Injecting Hierarchical Biological Priors into Graph Neural Networks for Flow
Cytometry Prediction

Fatemeh Nassajian Mojarrad  Lorenzo Bini  Thomas Matthes  Stéphane Marchand-
Maillet

###### Abstract

In the complex landscape of hematologic samples such as peripheral blood or
bone marrow derived from flow cytometry (FC) data, cell-level prediction
presents profound challenges. This work explores injecting hierarchical prior
knowledge into graph neural networks (GNNs) for single-cell multi-class
classification of tabular cellular data. By representing the data as graphs
and encoding hierarchical relationships between classes, we propose our
hierarchical plug-in method to be applied to several GNN models, namely, FCHC-
GNN, and effectively designed to capture neighborhood information crucial for
single-cell FC domain. Extensive experiments on our cohort of 19 distinct
patients, demonstrate that incorporating hierarchical biological constraints
boosts performance significantly across multiple metrics compared to baseline
GNNs without such priors. The proposed approach highlights the importance of
structured inductive biases for gaining improved generalization in complex
biological prediction tasks.

Machine Learning, ICML, Graph Neural Networks, Prior Knowledge, Biology

##  1 Introduction

Flow cytometry (FC) is a technology that provides rapid multi-parametric
analysis of single cells in solution. The expression of specific cell surface
and intracellular molecules is thereby detected by antibodies coupled to
fluorochromes. Exposing labelled cells to a laser leads to the emission of
fluorescence which is then collected by a detector. By analyzing more than
millions of cells in a short time, information can thereby be collected about
the presence or absence of up to fifty different cellular molecules per cell.
As a result, histologists obtain tabular data where each row corresponds to a
single cell and every column to a marker or feature (cell surface molecule –
Table 4, in Appendix A).

This technique is employed in medicine to assess the cellular composition of
complex body fluids like blood or bone marrow and to facilitate the diagnosis
of hematologic diseases such as leukemia. In collaboration with our
University’s Hospital, we have built a dataset obtained by the analysis of
bone marrow samples from 19 healthy/recovered patients with information
obtained up to $1^{\prime}512^{\prime}610$ cells/sample about the presence or
absence of twelve different surface molecules. Cellular processes are governed
by complex hierarchical relationships and neighborhood interactions that are
challenging to capture using traditional flat data representations. This
strongly motivates our decision to apply GNNs, which have emerged as a
powerful class of models capable of leveraging graph-structured data to reason
about relationships and dependencies. By representing biological entities and
their associations as nodes and edges in a graph, GNNs can effectively encode
neighborhood information and propagate signals across the graph topology
during training. However, even with this graph-based inductive bias, GNN
models may still struggle to generalize well without additional guided
constraints from biological domain knowledge.

In this work, we propose a novel approach to inject hierarchically-structured
priors (biological in our case) into the GNN framework for multi-class
prediction tasks on tabular data. Specifically, we encode the known
hierarchical relationships between different cell types or functional classes
as a tree-structured hierarchy imposed on the GNN output space. This
hierarchical constraint serves as an inductive bias that encourages the GNN to
respect the hierarchical dependencies inherent to the biological domain.
During training, we use a custom hierarchical loss function that accounts for
the hierarchical similarities between classes, in addition to the traditional
cross-entropy loss. This approach ensures that the model’s predictions not
only accurately classify instances into their respective leaf nodes (specific
cell types), but also respect the hierarchical groupings at higher levels of
the taxonomy (broader cell lineages or functional categories). Through
rigorous experiments on our real-world datasets, we demonstrate that
incorporating hierarchical biological priors significantly boosts the
performance of GNN models across a range of metrics. Our hierarchical GNN
approach achieves consistent gains over strong baselines that do not leverage
such structured domain knowledge. We provide in-depth analysis and empirical
evidence highlighting the importance of encoding inductive biases aligned with
the underlying data domains, especially for domains with rich hierarchical
relationships like cell-biology.

##  2 Related Work

The primary research focus in machine learning has predominantly been centered
on developing models for conventional classification problems, wherein an
object is assigned to a single class from a set of disjoint classes. However,
a distinct subset of tasks involves scenarios where classes are not disjoint
but rather organized hierarchically, giving rise to hierarchical
classification (HC). The hierarchical structure formalizing the interrelation
among classes can manifest as either a tree or a directed acyclic graph (DAG)
(Silla & Freitas, 2011). Hierarchical classification problems manifest across
a broad spectrum, spanning from musical genre classification (Ariyaratne &
Zhang, 2012; Iloga et al., 2018) to identification of COVID-19 in chest X-ray
images (Pereira et al., 2020), taxonomic classification of viral sequences in
metagenomic data (Shang & Sun, 2021) and text categorization (Javed et al.,
2021; Ma et al., 2022). Graph-based models have the capability to leverage
both the global and local characteristics inherent in networks relevant to
cell biology. GNNs have been recently adapted in various tasks, e.g. link
prediction (Zhang & Chen, 2018), graph classification (Duvenaud et al., 2015;
Lee et al., 2019), and node classification (Kipf & Welling, 2017; Velickovic
et al., 2017; Yang et al., 2016). In the context of HC problems, DiffPool
(Ying et al., 2018) propose a differentiable graph pooling module to generate
hierarchical representations of graphs using various GNNs in an end-to-end
fashion. Moreover, in the biological domains, HC-GNNs have been widely applied
for histological image classification (Hou et al., 2022), protein-protein
interaction (Gao et al., 2023), molecules’ properties prediction (Han et al.,
2023) and for cancer diagnosis and prognogis in digital pathology (Pati et
al., 2022).

##  3 Problem Description

In this research, we began by collecting raw data from the bone marrow of 19
patients who underwent a flow cytometric analysis for diagnostic purposes. All
samples were processed by the Diagnostics laboratory of the University
Hospital and no malignant disease was detected at this stage. Each patient’s
sample is primarily tabular data (see supplementary material for the futher
details) composed by a $N\times D$ matrix of values, where $D=12$ dimensions
(see Table 4 in Appendix A for details). It is therefore equivalent to a
$D$-dimensional point cloud that cytometrists visualize via 2D projections.
The result of such operations is the definition and quantification of
different cell populations according to their phenotype. Those populations of
interest were defined according to the markers expression, and then grouped
together into different categories as shown in Figure 1, serving as our main
biological hierarchical prior in this study.

for tree= edge=-¿, thick, l=1.5cm [Total cell population [CD45 pos
[Lymphocytes [B cells [Lambda pos] [Kappa pos] ] [NK cells] [T cells [CD8 T
cells] [CD4 T cells] ] ] [Monocytes] [Neutrophils] ] [CD45 neg] ]

Figure 1: Depiction of our HC set up.

The main focus is to classify different cell types present in the
heterogeneous samples into 12 classes, where 8 of them belong to leaf nodes
(green ones), while 4 of them are parental nodes (blue ones).

##  4 Archithecture of FCHC-GNN

Given a tabular dataset containing $n$ data samples and $m$ feature fields,
denoted by $X=[\mathbf{x}_{1},\cdots,\mathbf{x}_{n}]^{T}$, the $i$-th sample
in the table associated with $m$ feature values
$\mathbf{x}_{i}=[x^{(1)}_{i},\cdots,x^{(m)}_{i}]$ and a discrete label
$y_{i}$, part of a label hierarchy. Our learning goal is to find a mapping
function $f$ for which given $\mathbf{x}_{i}$, returns the predicted label
$\hat{y}_{i}$. In the rest of this section, we will describe our plug-in
method FCHC-GNN and the architecture of the network for the underlying
classification problem.

We consider a HC problem with a given set $\mathcal{C}$ of $C$ classes, which
are hierarchically organized as a directed tree (Figure 1 for example). If
there exists a path of non-negative length from class $A_{1}$ to class $A_{2}$
in the tree, we call $A_{2}$ a subclass of $A_{1}$. We assume to have a
mapping $f_{A}:\mathbb{R}^{m}\to[0,1]$ for every class $A$ such that
$\mathbf{x}\in\mathbb{R}^{m}$ is predicted to belong to $A$ whenever
$f_{A}(\mathbf{x})$ is greater than or equal to some threshold. When $A_{2}$
is a subclass of $A_{1}$, the model should guarantee that
$f_{A_{1}}(\mathbf{x})\geq f_{A_{2}}(\mathbf{x})$, for all
$\mathbf{x}\in\mathbb{R}^{m}$, in order to guarantee that the hierarchy
constraint is always satisfied independently from the threshold. Note that in
this model, every class is considered a subclass of itself. The case where
$f_{A_{1}}(\mathbf{x})<f_{A_{2}}(\mathbf{x})$ for some
$\mathbf{x}\in\mathbb{R}^{m}$ and $A_{2}$ is a subclass of $A_{1}$, is called
a hierarchy violation (Wehrmann et al., 2018). To simplify the notations, we
will write $f_{A}(\mathbf{x})$ as $f_{A}$, thus removing the functional
dependency of this term with respect to $\mathbf{x}$.

#### Overall structure

Given a set $\mathcal{C}=\\{A_{1},\cdots,A_{C}\\}$ of $C$ hierarchically
structured classes $A_{c}$, we first build a neural network composed of two
modules: a Early Module $\mathcal{H}$ with one output for each class in
$\mathcal{C}$ and an upper module, referred to as the Max Constraint Module
(MCM), consisting of a single layer that takes the output of the early module
as input and imposes the hierarchy constraint.

If $y_{A}$ is the ground truth label for class $A\in\mathcal{C}$ and the
prediction value $\mathcal{H}_{A}\in[0,1]$ is obtained by $\mathcal{H}$ for
class $A$, the MCM output for a class $A$ is given by:

$\text{MCM}_{A}=\max_{B\in\mathcal{S}_{A}}\mathcal{H}_{B}.$ (1)

where $\mathcal{S}_{A}$ is the set of subclasses of $A$ in $\mathcal{C}$. Note
that since the number of operations performed by $\text{MCM}_{A}$ is
independent from the depth of the hierarchy, for each class $A\in\mathcal{C}$,
the resulting model is scalable.

By construction, the following theorem therefore holds.

######  Theorem 4.1.

Let $\mathbf{x}\in\mathbb{R}^{m}$ be a data point. Let
$\mathcal{C}=\\{A_{1},\cdots,A_{C}\\}$ be a set of hierarchically structured
classes and let $\mathcal{H}$ be a early module with outputs
$\mathcal{H}_{A_{1}},\cdots,\mathcal{H}_{A_{C}}$
($\mathcal{H}_{A_{c}}\in[0,1]~{}\forall c$) given the input $\mathbf{x}$, and
$\text{MCM}_{A_{1}},\cdots,\text{MCM}_{A_{C}}$ be defined as in equation (1).
Then FCHC-GNNs do not admit hierarchy violations.

Second, in order to exploit the hierarchy constraint during training, FCHC-
GNNs are trained with max constraint loss (MCLoss), defined as follows

$\text{MCLoss}_{A}=-y_{A}\ln(\max_{B\in\mathcal{S}_{A}}y_{B}\mathcal{H}_{B})-(1-y_{A})\ln(1-\text{MCM}_{A}),$

for each class $A\in\mathcal{C}$. The global MCLoss is written as

$\text{MCLoss}=\sum_{A\in\mathcal{C}}\text{MCLoss}_{A}.$ (2)

The advantage of defining MCLoss instead of using the standard binary cross-
entropy loss is that the more ancestors a class has, the more likely it is
that FCHC-GNN trained with the standard binary cross-entropy loss would remain
stuck in spurious local optima. The MCLoss prevents this to happen.

As a result, FCHC-GNN has the ability of delegating the prediction for a class
$A_{i}$ to one of its subclasses $A_{j}$, thanks to MCM and the MCLoss. More
formally, the following proposition holds by construction.

######  Proposition 4.2.

Under the same assumptions as in Theorem 4.1, considering a class
$A_{i}\in\mathcal{C}$ and $A_{j}\in\mathcal{S}_{A_{i}}$ with $i\neq j$, FCHC-
GNN delegates the prediction on $A_{i}$ to $A_{j}$ for $\mathbf{x}$, if
$\text{MCM}_{A_{i}}=\mathcal{H}_{A_{j}}$ and
$\mathcal{H}_{A_{i}}<\mathcal{H}_{A_{j}}$.

#### Classification layer

We have designed the FCHC module to be compatible for GAT, SAGE and GCN, as
highlighted in our code publicly accessible at https://github.com111Anonymity
of the reviewing process., but due its general nature it can be easily
extended to broader GNNs. For sake of simplicity, in the following discussion
we explain how we built the Early Module $\mathcal{H}$ for the classification
as a Graph Attention Network (GAT), following (Velickovic et al., 2017). Given
$\\{\mathbf{x}_{i}\\}_{i=1}^{n}$ as nodes, and
$\mathbf{x}_{i}\in\mathbb{R}^{m}$, the node features are the input of the
graph attention layer. The output of the layer is a new set of node features
$\\{\mathcal{H}_{i}\\}_{i=1}^{n}$, where
$\mathcal{H}_{i}\in\mathbb{R}^{\tilde{m}}$, for which $\tilde{m}$ might be
different from $m$.

A GAT defines attention coefficients as importance scores of node
$\mathbf{x}_{j}$’s features to that of node $\mathbf{x}_{i}$. $W$ is the
weight matrix used to build attention coefficients. Masked attention
normalizes attention coefficients within every graph neighborhood
$\mathcal{N}(\mathbf{x}_{i})$. In our case, $\mathcal{N}(\mathbf{x}_{i})$
represents the Euclidean $k$-nearest neighbours of node $\mathbf{x}_{i}$. The
attention mechanism is a single-layer feedforward neural network,
parameterized by a weight vector $\vec{a}\in\mathbb{R}^{2\tilde{m}}$, over
which we apply the LeakyReLU non-linearity ( see Fig. 3(a) in Appendix C).

The normalized attention coefficients $\gamma_{ij}$ are then used to compute a
linear combination of the features corresponding to them, after applying a
nonlinearity $\sigma$, to serve as the final output features for each node.

We also move from single head attention to multi-head attention (with $L$
heads) by concatenating the output of all heads at every layer $i$ (but the
final layer). The expression for the final layer replaces concatenation by
averaging and can be written in our notation as

$\mathcal{H}_{i}=\sigma(\frac{1}{L}\sum_{l=1}^{L}\sum_{j\in\mathcal{N}(\mathbf{x}_{i})}\gamma_{ij}^{l}W^{l}\mathbf{x}_{j}).$
(3)

##  5 Experiments

In this section, we present the experimental results to illustrate the
behaviour of ours FCHC-GNNs for the multi-class classification problem on our
FC dataset. We also compare the proposed method with the respective flat
versions (meaning without the FCHC module attached), to highlight how the
injection of hierarchical prior knowledge really boost the performance for FC
classification with respect to flat structures. Moreover, since they are known
to be highly performant on tabular data, we also incuded Deep Neural Network
both with (FCHC-DNN) and without (DNN) hierarchical module. We have also
obtained competitive performance of our FCHC-GNN module on the public
Imagenet1000-(Mini) dataset, where we used the subcategory of ”animals”; these
additional results are presented in Appendix D due to space constraints.

Furthermore, in Appendix E as ablation studies, we further validate the
generality of our FCHC-GNN module by conducting experiments on other 30
different FC patients with a shallower depth of hierarchy. Across this cohort,
we observe that injecting the hierarchical priors becomes necessary as the
hierarchy grows deeper and more complex. For shallow hierarchies, the
performance gains are modest, but they become substantially larger when
dealing with intricate, multi-level hierarchies where using prior knowledge
becomes fundamental.

Due to space constraints, see Appendix A for experimental setup, metrics
definition and graph-construction details. Please see also Appendix E for
model’s ablation studies on different hierarchy scenario and hyperparameters
details. Our hierarchical approach achieves consistent gains over same
architectures that do not leverage such structured domain knowledge, as shown
in Table 1, 2 and 3, where we compute hierarchical-metrics and the correct
predicted classes across all the patients in the dataset. All the results are
within the $\text{std}=\pm 0.1$, and experiments have been averages across
four different seeds.

Table 1: Average metrics across all patients using different models for FCHC
module.
Metrics | FCHC-GAT | FCHC-SAGE | FCHC-GCN | FCHC-DNN
---|---|---|---|---
hp | 0.88 | 0.92 | 0.83 | 0.74
hr | 0.85 | 0.83 | 0.81 | 0.83
hf | 0.87 | 0.88 | 0.82 | 0.78
Table 2: Average ratios of correct predicted classes across all patients using
different models for FCHC module.
Label | FCHC-GAT | FCHC-SAGE | FCHC-GCN | FCHC-DNN
---|---|---|---|---
CD45 pos | 100.00 | 100.00 | 100.00 | 81.12
CD45 neg | 100.00 | 100.00 | 100.00 | 83.43
Lymphocytes | 85.21 | 97.91 | 82.14 | 81.06
Monocytes | 90.31 | 93.42 | 95.52 | 90.07
Neutrophils | 80.15 | 93.38 | 75.26 | 80.84
B Cells | 94.86 | 97.49 | 98.52 | 85.83
NK Cells | 97.48 | 97.50 | 97.98 | 94.01
T Cells | 83.00 | 97.84 | 86.74 | 92.07
Lambda pos | 97.77 | 98.88 | 96.52 | 86.45
Kappa Pos | 96.56 | 98.61 | 99.21 | 98.00
CD8 T Cells | 90.15 | 95.26 | 94.25 | 91.99
CD4 T Cells | 89.81 | 91.63 | 92.07 | 90.06

###  5.1 Discussion of Our Results

The results presented in this paper underscore the novelty and effectiveness
of our newly introduced FCHC-GNN framework, in performing node classification
tasks on FC dataset. However, the presented approach is rather general and can
be easily extend to broader applications of GNN, in case there is the need to
leverages the inherent hierarchical structure of the data, a feature that
frequently appears in real-word dataset. Indeed, as shown in Table 3 where the
FCHC model has not been applied, injecting hierarchical biological priors
allows the models to not overlook any of the leaf nodes during the prediction.
As an example, in Table 2 and 3 we clearly see how the performances drops from
$91.63\%$ with hierarchical priors to a $67.80\%$ without, on CD4 T cells for
GraphSAGE model (Hamilton et al., 2017). Similarly, for NK cells the correct
predicted labels drammatically goes down from $97.48\%$ to $0.15\%$ for the
GAT, highlighting once again the need of including the FCHC module to caputure
the biology complexity of these cells population.

Table 3: Average ratios of corrected predicted classes across all patients
without using the FCHC module.
Label | GAT | SAGE | GCN | DNN
---|---|---|---|---
CD45 neg | 34.94 | 98.25 | 97.10 | 98.44
Monocytes | 0.75 | 1.46 | 0.50 | -
Neutrophils | 100.00 | 98.96 | 99.73 | 98.77
NK Cells | 0.15 | 1.57 | - | -
Lambda pos | - | 0.88 | 0.07 | 1.40
Kappa Pos | 0.24 | - | 0.10 | 0.70
CD8 T Cells | - | 20.41 | 13.68 | 1.01
CD4 T Cells | 0.20 | 67.80 | 40.08 | 31.95

Given the classification power of FCHC module, we took a deeper look at each
model’s performance, and to allow better visualization of the predicted cell
patterns, we picked one random patient and we plotted its t-SNE embeddings
(Van der Maaten & Hinton, 2008), as shown in Figure 5 of Appendix G. Moreover,
we further investigated the interpretability of our module through the
expression of feature importance, as resulted in Appendix B.

##  6 Conclusion and Future Work

We have presented a novel hierarchical GNN framework, namely FCHC-GNN, that
injects structured biological priors to enhance multi-class prediction on FC
cellular data. By encoding known hierarchical relationships between cell types
and functional classes, our model effectively captures the rich dependencies
innate to biological domains while operating on a tabular input
representation. Extensive experiments across ours FC dataset demonstrate
substantially improved performance over strong baselines that fail to leverage
such hierarchical constraints. Moreover, as demonstrated over well-know
ImageNet data, our FCHC-GNN module is easily extendible to other domain where
classes may be organized within an hierarchy.

The consistent gains validate our hypothesis that tailored inductive biases
aligned with the underlying tabular data distributions are crucial for
achieving better generalization. FCHC-GNN architectures introduces a novel
method to enhance learning on tabular inputs with prior domain knowledge
expressed as hierarchies. Looking ahead, this general approach holds promise
for biomedical learning tasks beyond single cell classification.

## Acknowledgements

This work is partly funded by the Swiss National Science Foundation under
grant number 207509 ”Structural Intrinsic Dimensionality”.

## References

  * Ariyaratne & Zhang (2012) Ariyaratne, H. B. and Zhang, D.  A novel automatic hierachical approach to music genre classification.  In _Proceedings of the IEEE International Conference on Multimedia and Expo Workshops_ , pp. 564–569. IEEE, July 2012.
  * Duvenaud et al. (2015) Duvenaud, D. K., Maclaurin, D., Iparraguirre, J., Bombarell, R., Hirzel, T., Aspuru-Guzik, A., and Adams, R. P.  Convolutional networks on graphs for learning molecular fingerprints.  _Advances in neural information processing systems_ , 28, 2015.
  * Gao et al. (2023) Gao, Z., Jiang, C., Zhang, J., Jiang, X., Li, L., Zhao, P., Yang, H., Huang, Y., and Li, J.  Hierarchical graph learning for protein–protein interaction.  _Nature Communications_ , 14(1):1093, 2023.
  * Hamilton et al. (2017) Hamilton, W., Ying, Z., and Leskovec, J.  Inductive representation learning on large graphs.  _Advances in neural information processing systems_ , 30, 2017.
  * Han et al. (2023) Han, S., Fu, H., Wu, Y., Zhao, G., Song, Z., Huang, F., Zhang, Z., Liu, S., and Zhang, W.  Himgnn: a novel hierarchical molecular graph representation learning framework for property prediction.  _Briefings in Bioinformatics_ , 24(5):bbad305, 2023.
  * Hou et al. (2022) Hou, W., Huang, H., Peng, Q., Yu, R., Yu, L., and Wang, L.  Spatial-hierarchical graph neural network with dynamic structure learning for histological image classification.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pp. 181–191. Springer, 2022.
  * Iloga et al. (2018) Iloga, S., Romain, O., and Tchuenté, M.  A sequential pattern mining approach to design taxonomies for hierarchical music genre recognition.  _Pattern Analysis and Applications_ , 21(2):363–380, May 2018.
  * Javed et al. (2021) Javed, T. A., Shahzad, W., and Arshad, U.  Hierarchical text classification of urdu news using deep neural network.  https://arxiv.org/abs/2107.03141, 2021.
  * Kipf & Welling (2017) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  https://openreview.net/pdf?id=SJU4ayYgl, 2017.
  * Kiritchenko et al. (2006) Kiritchenko, S., Matwin, S., Nock, R., and Famili, A. F.  Learning and evaluation in the presence of class hierarchies: Application to text categorization.  In _Proceedings of the Conference of the Canadian Society for Computational Studies of Intelligence_ , pp. 395–406, Québec, Canada, June 2006. Springer Berlin Heidelberg.
  * Lee et al. (2019) Lee, J., Lee, I., and Kang, J.  Self-attention graph pooling.  In _International conference on machine learning_ , pp. 3734–3743. PMLR, 2019.
  * Ma et al. (2022) Ma, Y., Liu, X., Zhao, L., Liang, Y., Zhang, P., and Jin, B.  Hybrid embedding-based text representation for hierarchical multi-label text classification.  _Expert Systems with Applications_ , 187:115905, May 2022.
  * Pati et al. (2022) Pati, P., Jaume, G., Foncubierta-Rodriguez, A., Feroce, F., Anniciello, A. M., Scognamiglio, G., Brancati, N., Fiche, M., Dubruc, E., Riccio, D., et al.  Hierarchical graph representations in digital pathology.  _Medical image analysis_ , 75:102264, 2022.
  * Pereira et al. (2020) Pereira, R. M., Bertolini, D., Teixeira, L. O., Jr, C. N. S., and Costa, Y. M.  Covid-19 identification in chest x-ray images on at and hierarchical classification scenarios.  _Computer Methods and Programs in Biomedicine_ , 194:105532, October 2020.
  * Shang & Sun (2021) Shang, J. and Sun, Y.  Cheer: hierarchical taxonomic classification for viral metagenomic data via deep learning.  _Methods_ , 189:95–103, May 2021.
  * Silla & Freitas (2011) Silla, C. N. and Freitas, A. A.  A survey of hierarchical classifcation across different application domains.  _Data Mining and Knowledge Discovery_ , 22(1):31–72, 2011.
  * Van der Maaten & Hinton (2008) Van der Maaten, L. and Hinton, G.  Visualizing data using t-sne.  _Journal of machine learning research_ , 9(11), 2008.
  * Velickovic et al. (2017) Velickovic, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., and Bengio, Y.  Graph attention networks.  _stat_ , 1050(20):10–48550, 2017.
  * Wehrmann et al. (2018) Wehrmann, J., Cerri, R., and Barros, R.  Hierarchical multi-label classification networks.  In _International conference on machine learning_ , pp. 5075–5084. PMLR, July 2018.
  * Yang et al. (2016) Yang, Z., Cohen, W., and Salakhudinov, R.  Revisiting semi-supervised learning with graph embeddings.  In _International conference on machine learning_ , pp. 40–48. PMLR, June 2016.
  * Ying et al. (2018) Ying, Z., You, J., Morris, C., Ren, X., Hamilton, W., and Leskovec, J.  Hierarchical graph representation learning with differentiable pooling.  _Advances in neural information processing systems_ , 31, 2018.
  * Zhang & Chen (2018) Zhang, M. and Chen, Y.  Link prediction based on graph neural networks.  _Advances in neural information processing systems_ , 31, 2018.

##  Appendix A Description of Our Dataset

The raw data obtained from the cytometer are saved in Flow Cytometry Standard
(FCS) files. Each FCS file comprises multidimensional data that corresponds to
more than one million of cells, with each cell characterized by various
parameters such as size, granularity, and fluorescence intensity, giving them
the structure of tabular data (see Table 4 for a complete list of markers in
our data).

Table 4: Flow cytometry data markers.
Marker # | Marker | Explanation
---|---|---
0 | FS INT | Forward Scatter (FSC) - Cell’s size
1 | SS INT | Side Scatter (SSC) - Cell’s granularity
2 | CD14-FITC | Cluster of Differentiation 14 - Antigen
3 | CD19-PE | Cluster of Differentiation 19 - Antigen
4 | CD13-ECD | Cluster of Differentiation 13 - Antigen
5 | CD33-PC5.5 | Cluster of Differentiation 33 - Antigen
6 | CD34-PC7 | Cluster of Differentiation 34 - Antigen
7 | CD117-APC | Cluster of Differentiation 117 - Antigen
8 | CD7-APC700 | Cluster of Differentiation 7 - Antigen
9 | CD16-APC750 | Cluster of Differentiation 16 - Antigen
10 | HLA-PB | Human Leukocyte Antigen
11 | CD45-KO | Cluster of Differentiation 45 - Antigen

The synthesis of machine learning and biological domain knowledge holds
potential for advancing scientific understanding and biomedical applications.
While machine learning models excel at finding intricate patterns in data,
their performance can be limited without the inductive biases, in particular
for tabular data which typically lacks of it, and constraints derived from
prior biological principles.

###  A.1 Hierarchical Metrics Definition and Graph-Construction

In our experiments, we construct the $k$-nearest neighbors graphs with $k=7$
neighbors. We perform a 7-fold test procedure where 1/7th of the data (2 or 3
patients out of 19) is held out as a test set and the rest (6/7th of data) as
a training set. We then perform a 4-fold cross validation procedure within the
training set to set the hyperparameters.

The metrics are reported using the extensions of the renowned metrics of
precision, recall and F-score, but tailored to the HC setup. We implement the
metrics of hierarchical precision ($hp$), hierarchical recall ($hr$) and
hierarchical F-score ($hf$) defined by (Kiritchenko et al., 2006):

$hp=\frac{\sum_{i}|\alpha_{i}\cap\beta_{i}|}{\sum_{i}|\alpha_{i}|},~{}hr=\frac{\sum_{i}|\alpha_{i}\cap\beta_{i}|}{\sum_{i}|\beta_{i}|},~{}hf=\frac{2\times
hp\times hr}{hp+hr}.$

(4)

where $\alpha_{i}$ is the set consisting of the most specific classes
predicted for test sample $i$ and all their ancestor classes and $\beta_{i}$
is the set containing the true most specific classes of test sample $i$ and
all their ancestors.

##  Appendix B Interpretability and Feature Importance Analysis

An important aspect of our FCHC plug-in module is its interpretability, as
demonstrated in Figure 2 for the FCHC-GAT version.

Figure 2: Feature importance for FCHC-GAT attributed by $\mathcal{H}$. Feature
labels correspond to that of Table 4.

This feature importance analysis not only provides insights into the decision-
making process of our model but also paves the way for potential improvements
and adaptations of our model for other applications. As we can see, the
features that are given the most importance (at least above 40%) are 2,1,5 and
0, therefore we broke down the choice of the model for each one. In order of
importance, we had:

  * •

CD14-FITC: This feature reflects the level of CD14 expression, a receptor
found on the cell surface that interacts with lipopolysaccharide (LPS), a
component of bacterial cell walls. Predominantly found on monocytes,
macrophages, and activated granulocytes, CD14 plays a pivotal role in
mediating the immune response to bacterial infections. The high weightage
given to this feature suggests that our model is adept at distinguishing
between cells involved in innate immunity and those that aren’t, thereby
potentially detecting the presence of bacterial infection in the sample.

  * •

Side Scatter (SSC)-Cell’s granularity: This feature measures the level of
light scatter at a 90-degree angle to the laser beam, reflecting the internal
complexity or granularity of the cell. This could include the presence of
granules, nuclei, or other organelles. The high weightage of this feature
implies that our model can differentiate between cells of varying complexity,
such as lymphocytes, monocytes, and granulocytes, or identify cells with
abnormal granularity, such as blast cells or malignant cells.

  * •

CD33-PC5.5: This feature indicates the expression level of CD33, a cell
surface receptor of the sialic acid-binding immunoglobulin-like lectin
(Siglec) family. CD33 is expressed by myeloid cells like monocytes,
macrophages, granulocytes, and mast cells, and it modulates the immune
response by inhibiting the activation of these cells. The high weightage of
this feature suggests that our model can distinguish between myeloid and non-
myeloid cells, or detect the expression of CD33 as a marker for certain types
of leukemia, such as acute myeloid leukemia (AML) or chronic myelomonocytic
leukemia (CMML).

  * •

Forward scatter (FSC)-Cell’s size: This feature measures the level of light
scatter along the path of the laser beam, which is proportional to the
diameter or surface area of the cell. This can be used to discriminate cells
by size. The importance placed on this feature suggests that our model can
differentiate between cells of different sizes, such as small lymphocytes and
large monocytes.

##  Appendix C Graphical Representation of Multi-Head Attention

Figure 3 works as recap of the structural representation for the applied
multi-head attention.

Figure 3: (a): Computing the normalized attention coefficients $\gamma_{ij}$
(b): Multi-head attention of node 1 on its neighborhood. Arrows show
concatenation or averaging of attention (adapted from (Velickovic et al.,
2017)).

##  Appendix D ImageNet1000-Mini Results

Table 5 shows competitive performance of our FCHC-GNN module on hierarchical
image classification task as well. For sake of simplicity, we only showed our
plug-in module being applied to a GAT, against its normal version and the
average state-of-art performances (accordingly to what written on the
benchmark https://paperswithcode.com/sota/image-classification-on-imagenet).

Table 5: Average accuracy percentage of SOTA methods on Imagenet1000-(Mini).
Results are averaged with a $\text{std}=\pm 0.01$ across four different seeds.
Model | Imagenet1000-Mini
---|---
FCHC-GAT (Ours) | 91.88
GAT | 3.91
Average SOTA |  88.50 - 92.40

##  Appendix E Ablation Studies on a Shallow Experiments Setup

To compare the results given in Table 1 \- 3, we present the results on the
second cohort of FC patients, provided by our Univerisity Hospitaly, where 30
different bone-marrows samples have been collected and anylized in the same
way as the previous nineteens. The depth of hierarchy is shallower than that
of the main one, where only one step of parental nodes have been considered,
as depicted in Figure 4. As shown in Table 7 and Table 9 on shallow
hierarchies, the problem more closely resembles a flat classification task,
and thus baseline GNN performance remains relatively good. This is especially
true for datasets with large populations of certain cell types, where the
inherent class imbalance acts as an inductive bias, as depicted from the
hierarchical F-Score results on Table 8. However, when the hierarchies grow
deeper and more complex, the use of our FCHC-GNN module becomes strictly
necessary to obtain better interpretability and generalization. For example,
Table 6 already shows a consistant improvement on the hF-Score metrics with
respect to standard one without the hierarchical module, mainly due to the
correct prediction of Myeloid HSPC and other minor HSPC cells.

Deeper hierarchies exacerbate the issue of conflicting gradients between nodes
at different hierarchy levels during training. Without explicit hierarchical
constraints, as in the MCLoss, GNNs struggles to simultaneously optimize for
accurate predictions at the leaf node level while respecting parental
groupings (higher-level). Our hierarchical module resolves this tension
without suffering of extra time-consumation problem as explained in Appendix
F.3, allowing attached models to leverage hierarchical priors as a strong
inductive bias for learning more interpretable representations aligned with
the known biological taxonomy.

for tree= edge=-¿, thick, l=1.5cm [Root [T cells] [B cells] [Monocytes] [Mast
cells] [HSPC [Myeloid HSPC] [Lymphoid HSPC] [Other HSPC] ] ]

Figure 4: Depiction of the shallow FC setup. Table 6: Average metrics across
30 patients using different models for FCHC module.
Metrics | FCHC-GAT | FCHC-SAGE | FCHC-GCN | FCHC-DNN
---|---|---|---|---
hp | 0.94 | 0.97 | 0.95 | 0.87
hr | 0.93 | 0.97 | 0.94 | 0.83
hf | 0.94 | 0.97 | 0.94 | 0.85
Table 7: Average ratios of correct predicted classes across 30 patients using
different models for FCHC module.
Label | FCHC-GAT | FCHC-SAGE | FCHC-GCN | FCHC-DNN
---|---|---|---|---
T cells | 93.87 | 98.85 | 94.49 | 75.04
B cells | 95.45 | 98.28 | 96.03 | 86.82
Monocytes | 98.12 | 99.21 | 98.57 | 95.24
Mast cells | 99.79 | 99.93 | 99.79 | 99.79
HSPC | 95.03 | 97.85 | 95.73 | 91.03
Myeloid HSPC | 94.54 | 96.51 | 94.79 | 93.40
Lymphoid HSPC | 97.40 | 98.17 | 97.40 | 97.40
Other HSPC | 99.46 | 99.46 | 99.46 | 99.46
Table 8: Average metrics across 30 patients without using the FCHC module.
Metrics | GAT | SAGE | GCN | DNN
---|---|---|---|---
Precision | 0.92 | 0.84 | 0.92 | 0.84
Recall | 0.76 | 0.74 | 0.81 | 0.71
F1 Score | 0.75 | 0.64 | 0.79 | 0.61
Table 9: Average ratios of corrected predicted classes across 30 patients
without using the FCHC module.
Label | GAT | SAGE | GCN | DNN
---|---|---|---|---
T cells | 81.83 | 98.72 | 90.14 | 96.91
B cells | 51.50 | 10.96 | 57.22 | 9.87
Monocytes | 82.89 | 69.26 | 80.89 | 56.19
Mast cells | 98.22 | 86.51 | 52.62 | 53.33
Myeloid HSPC | 26.16 | - | 12.37 | -
Lymphoid HSPC | 90.75 | 12.96 | 84.50 | 26.43
Other HSPC | 23.11 | - | 32.48 | -

##  Appendix F FCHC Plug-In Module Implementation

This section provides additional details on how the hierarchical constraint
module was implemented in all the tested models. For sake of simplicity, we
choose to illustrate how we designed the FCHC-GAT module, while for other
FCHC-GNNs implementations refer to our GithHub repository https://github.com.
Here are the key steps:

  * •

It first defines the hierarchical structure by creating a directed graph $g$
from the class hierarchy string ATTRIBUTE_class, as shown in Listing
LABEL:code:MCM below. This maps out the parent-child relationships between
classes.

  * •

It converts this graph to a matrix $R$ that encodes the ancestor-descendant
relationships, where $R[i,j]=1$ if class $i$ is descendant of class $j$, $0$
otherwise.

  * •

Then the model includes GAT layers, where each of them does message passing
and computes attention coefficients between neighboring nodes.

  * •

During training, the raw GAT predictions are returned.

  * •

During inference, the predictions are passed through a constraint function
get_constr_out that takes the max over descendants to enforce the hierarchy,
as depicted in Listing LABEL:code:leukograph below.

  * •

The loss function is a modified cross-entropy loss MCLoss that also lets
subclasses inform about their parent classes predictions.

  * •

The model takes a graph and node features as input, does message passing and
attention in the GAT layers, applies the hierarchical constraint, and outputs
a prediction vector for each node.

In summary, it implements a hierarchical classifier using a GAT augmented with
constraints and loss terms to incorporate the class hierarchy information. The
hierarchy structure is predefined and encoded in the constraint matrix $R$.

###  F.1 Constraint Layer

The hierarchy consistency is imposed through a constraint layer added on top
of the base classifier network. This constraint layer, referred to as the Max
Constraint Module (MCM), takes the output predictions from the base network
and ensures the hierarchical constraint is satisfied. The MCM output for any
class $A$ is computed as:

$\displaystyle\text{MCM}_{A}=\max(\mathcal{H}_{B}\text{ where }B\text{ is
subclass of }A),$ (5)

where $\mathcal{H}_{B}$ is the prediction value from the base network for
class $B$. This takes the maximum over all subclasses of $A$, guaranteeing
that the prediction for class $A$ will be greater than or equal to the
prediction for any of its subclasses. By constructing the output this way, the
hierarchy constraint is inherently satisfied regardless of the threshold used
downstream.

⬇

1# List all the classes

2to_skip = [’root’]

3ATTRIBUTE_class =
"1,1_1,1_1_1,1_1_1_1,1_1_1_2,1_1_2,1_1_3,1_1_3_1,1_1_3_2,1_2,1_3,2"

4

5# Store nodes and direct connections

6g = nx.DiGraph()

7for branch in ATTRIBUTE_class.split(’,’):

8 term = branch.split(’_’)

9 if len(term)==1:

10 g.add_edge(term[0], ’root’)

11 else:

12 for i in range(2, len(term) + 1):

13 g.add_edge(’.’.join(term[:i]), ’.’.join(term[:i-1]))

14

15nodes = sorted(g.nodes(), key=lambda x: (len(x.split(’.’)),x))

16

17AA = np.array(nx.to_numpy_array(g, nodelist=nodes))

18

19# Compute matrix of ancestors R

20# Given C classes, R is an (C x C) matrix

21# where R_ij = 1 if class i is descendant of class j

22R = np.zeros(AA.shape)

23np.fill_diagonal(R, 1)

24gg = nx.DiGraph(AA) # AA is the matrix where the direct connections are
stored

25for i in range(len(AA)):

26 ancestors = list(nx.descendants(gg, i)) # Need to use the function
nx.descendants()

27 #because in the directed graph the edges have source

28 #from the descendant and point towards the ancestor

29 if ancestors:

30 R[i, ancestors] = 1

31R = torch.tensor(R)

32#Transpose to get the descendants for each node

33R = R.transpose(1, 0)

34R = R.unsqueeze(0).to(device)

Listing 1: Python code snippet generating a matrix of ancestors $R$ based on
the class relationships in the directed graph. Each class is represented as a
node, and edges define parent-child connections.

###  F.2 Constrained Loss Function

The loss function used during training is the max constraint loss (MCLoss),
defined as:

$\displaystyle\text{MCLoss}_{A}=-y_{A}\log\left(\max_{B\text{ in
subclasses\\_of\\_A}}(y_{B}\cdot\mathcal{H}_{B})\right)-(1-y_{A})\log(1-\text{MCM}_{A}),$
(6)

where $y_{A}$ is the ground truth label. This loss allows the model to
leverage the hierarchical information by letting predictions from subclasses
inform the training for parent classes.

###  F.3 Time Complexity

The time complexity of the constraint layer is $O(C)$ where $C$ is the number
of classes, since it just takes a maximum over the subclasses for each class.
So it scales linearly with the size of the hierarchy and does not depend on
the depth. The overall model maintains the same asymptotic time complexity as
the base classification network used, with just an additional linear factor
for the constraint layer. Hence, the hierarchical structure does not add
significant overhead during training or inference.

⬇

1def get_constr_out(x, R):

2 ’’’ Given the output of the graph neural network x returns the output of

3 MCM given the hierarchy constraint expressed in the matrix R ’’’

4 c_out = x.double()

5 c_out = c_out.unsqueeze(1)

6 c_out = c_out.expand(len(x),R.shape[1], R.shape[1])

7 R_batch = R.expand(len(x),R.shape[1], R.shape[1])

8 final_out, _ = torch.max(R_batch*c_out.double(), dim = 2)

9 return final_out

10

11# Define the FCHC plug-in module

12class FCHC_GAT(nn.Module):

13 ’’’ During training it returns the not-constrained output that is then
passed to MCLoss’’’

14 def __init__(self,R):

15 super(FCHC_GAT, self).__init__()

16 self.R = R

17 self.nb_layers = 2

18 self.num_heads = 2

19 self.out_head = 2

20 self.hidden_dim = 32

21 self.input_dim = 12

22 self.output_dim= len(set(ATTRIBUTE_class.split(’,’)))+1 # We do not
evaluate

23 #the performance of the model on the ’root’ node

24 gat_layers = []

25 for i in range(self.nb_layers):

26

27 if i == 0:

28 gat_layers.append(GATLayer(self.input_dim, self.hidden_dim,

29 self.num_heads,True,0.4))

30 elif i == self.nb_layers - 1:

31 gat_layers.append(GATLayer(self.hidden_dim*self.num_heads,

32 self.output_dim, self.out_head, False,0.2))

33 else:

34 gat_layers.append(GATLayer(self.hidden_dim*self.num_heads,

35 self.hidden_dim, self.num_heads, True,0.2))

36 self.gat_layers = nn.ModuleList(gat_layers)

37

38 self.sigmoid = nn.Sigmoid()

39 self.f = nn.ReLU()

40 self.reset_parameters()

41 self.drop = nn.Dropout(0.2)

42

43 def reset_parameters(self):

44 for gat_layer in self.gat_layers:

45 gat_layer.reset_parameters()

46

47 def forward(self, data):

48 x, edge_index= data.x, data.edge_index

49 for i in range(self.nb_layers):

50 x = self.gat_layers[i](x, edge_index)

51 if i != self.nb_layers - 1:

52 x = self.f(x)

53 x = self.drop(x)

54 else:

55 x= self.sigmoid(x)

56 if self.training:

57 constrained_out = x

58 else:

59 constrained_out = get_constr_out(x, self.R )

60 return constrained_out

Listing 2: FCHC-GAT plug-in module definition with hierarchical constraints.
The model incorporates the hierarchy matrix $R$ to enforce class relationships
during training. After performing message passing and attention computations,
the model’s output is constrained during inference using the get_constr_out
function.

##  Appendix G t-SNE Embedding Visualization of FCHC-GAT Module

An important aspect of our FCHC-GAT model that warrants special attention is
its ability to accurately capture and classify even small populations of
certain cell classes. This is a significant achievement, as demonstrates the
recall and precision of our model. The ability of FCHC-GAT is further
illustrated through visualization of a particular sub-population of cells,
particularly interested for physicians. As shown in Figure 5, the t-SNE
scatter plots provides a visual testament to the model’s ability to accurately
classify different cell types of interested, in this case NK Cells, CD4 T
Cells, Kappa Pos, Lambda Pos, Monocytes, Neutrophils and CD8 T cells.

Figure 5: Embeddings t-SNE projection for patient 7, for the FCHC-GAT module.