# Signature of the Maximally Symmetric 2HDM via $W^{\pm}/Z$-Quadruplet
Productions at the LHC

N. Darvishi∗,† and M.R. Masouminia‡ ∗School of Physics and Astronomy, The
University of Manchester, Manchester M13 9PL, United Kingdom
†Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland
‡Institute for Particle Physics Phenomenology, Department of Physics, Durham
University, Durham DH1 3LE, United Kingdom

###### Abstract

ABSTRACT

We consider the Maximally Symmetric Two-Higgs Doublet Model (MS-2HDM) in which
the so-called Standard Model (SM) alignment can be achieved naturally by the
virtue of an SO(5) symmetry imposed on the 2HDM. We investigate the signature
of the MS-2HDM via $pp\to HX\to VV^{*}X$ and $pp\to HHX\to
VV^{*}V^{\prime}V^{\prime*}X$ processes at the LHC for different values of
$\tan\beta$. We perform our calculations with NLO QCD accuracy, using the
Herwig 7 multi-purpose event generator at $\sqrt{s}=13$ TeV center-of-mass
energy. We show that the production of single SM-like Higgs bosons via
$W^{\pm}/Z$-pairs is completely aligned with the SM. Interestingly, the
presence of the heavy Higgs states significantly enhances the cross-section
for the $W^{\pm}/Z$-quadruplet production channels in the low-$p_{\perp}$
regions. These vital analyses may aid the future discovery of this minimal and
very predictive extension of the SM and can be generalised to other
realisations of the 2HDM.

##  I Introduction

For many years, the Standard Model (SM) has been the cornerstone for our
understanding of the fundamental interactions of Particle Physics
Glashow:1961tr ; Goldstone:1962es ; Weinberg:1967tq ; Salam:1968 . This was
brought to its climax with the discovery of the Higgs boson at CERN’s Large
Hadron Collider (LHC) Englert:1964et ; Higgs:1964pj . The data collected from
this discovery imposes constraints over the coupling strengths of the Higgs
boson, primarily to the electroweak (EW) gauge bosons ($V=W^{\pm},\;Z$), which
are very close to SM predictions ATLAS:2019aqa ; CMS:2019vxb . Despite all
these achievements, the SM still falls short of answering some of the most
profound questions such as the origin of the observed matter-antimatter
asymmetry and the dark matter relic abundance in the Universe. This has fueled
numerous theoretical and experimental scrutinise in the study of theories
Beyond the SM (BSM), particularly for models with extended Higgs sectors.
These new-born theories regardless of their motivational and structural
differences must restore those predictions of the SM that are consistent with
the LHC observations. This is only possible within the so-called SM alignment
limit Ginzburg:1999fb ; Chankowski:2000an ; Delgado:2013zfa ; Carena:2013ooa ;
Dev:2014yca ; Bernon:2015qea ; Darvishi:2016tni ; Darvishi:2017bhf ;
Benakli:2018vqz ; Lane:2018ycs .

One of the simplest extensions of the SM is the Two-Higgs Doublet Model
(2HDM), which enriches the SM scalar sector by introducing a second complex
scalar doublet Lee:1973iz ; Pilaftsis:1999qt ; Ginzburg:1999fb ; Branco:2005em
; Delgado:2013zfa . This extension can provide new sources of CP violation
Lee:1973iz ; Pilaftsis:1999qt , introduce stable scalar DM candidates
Silveira:1985rk ; Bonilla:2014xba ; Krawczyk:2015xhl , and give rise to EW
baryogenesis Kuzmin:1985mm ; Cohen:1993nk . Interestingly, the potential of
this model contains the maximum number of $13$ distinct SU(2)L-preserving
accidental symmetries as subgroups of the maximal symmetry Sp(4) $\sim$ SO(5)
Pilaftsis:2011ed ; Battye:2011jj . Thereby, the most minimal version of the
2HDM is an SO(5)-invariant potential, the so-called Maximally Symmetric 2HDM
(MS-2HDM) Pilaftsis:2011ed ; Battye:2011jj ; Dev:2014yca ; Pilaftsis:2016erj ;
Dev:2017org ; Hanson:2018uhf ; Darvishi:2019ltl ; Darvishi:2019dbh ;
Darvishi:2020teg ; Birch-Sykes:2020btk . In MS-2HDM, the SM alignment can
emerge naturally as a consequence of an accidental SO(5) symmetry in the Higgs
sector. This symmetry can be broken explicitly by the renormalization group
(RG) effects and softly by the bilinear scalar mass term $m^{2}_{12}$. A
remarkable feature of this model is that all quartic couplings can unify at
very large scales $\mu_{X}\sim 10^{9}\,$–$\,10^{20}$ GeV, for a wide range of
$\tan\beta$ values and charged Higgs-boson masses Darvishi:2019ltl ;
Darvishi:2020teg . This unique feature aids to gain a minimal and very
predictive model which is governed only by three parameters: the quartic
coupling unification scale $\mu_{X}$, the mass of charged Higgs $M_{h^{\pm}}$
(or $m^{2}_{12}$) and the value of $\tan\beta$. These three parameters allow
one to determine the entire Higgs-mass spectrum of the model.

The production of the EW gauge vector boson pairs and quadruplets have always
been among the critical observations in the on-going attempts to probe the SM
and search for signs of BSM physics at the LHC, within its high-energy
hadronic scattering data. Also, these events have played a key role in the LHC
precision measurements as well as the estimation of the irreducible
backgrounds in Higgs boson searches. Moreover, since these processes are
amongst the largest Higgs-tagged signatures at the current LHC energies,
observing a distinct deviation from the SM theoretical predictions may be
directly interpreted as a signal for BSM physics Darvishi:2016fwo ;
Darvishi:2019uzp . Therefore, it would be interesting to look for the possible
signature of the MS-2HDM at the LHC via $W^{\pm}/Z$-pair and -quadruplet
production events. In this paper, firstly we ensure that the predictions of
the MS-2HDM for the $pp\to HX\to VV^{*}X$ production rates are aligned with
their SM counterparts. Then, we evaluate the signature of the MS-2HDM via
$pp\to HHX\to VV^{*}V^{\prime}V^{\prime*}X$ events at the LHC for different
values of $\tan\beta$. We calculate these production rates up to one QCD loop
and 2 jets, using Herwig 7 (v7.2.1) event generator Bahr:2008pv ;
Bellm:2015jjp ; Bellm:2017bvx ; Bellm:2019zci . The corresponding amplitudes
are provided by MadGraph5 (v2.7.3) Alwall:2014hca  and matched to the NLO
corrections using Matchbox Platzer:2011bc ; Bellm:2019wrh . The produced
underlying events are showered by an AO MC@NLO matched QCD+QED+EW111The
QCD+QED+EW parton shower scheme is a new addition to Herwig 7 which is
introduced in Masouminia:2021  and will be available to the public with the
v7.3.0 release. parton shower Frixione:2002ik ; Masouminia:2021 . Finally, the
results of these simulations have been analyzed using Rivet (v3.1.1)
Buckley:2010ar .

The layout of the paper is as follows. After this introductory section,
Section II briefly reviews the basic features of the 2HDM and the naturally
aligned MS-2HDM. We also outline the Higgs-mass spectrum and our misalignment
predictions for Higgs-boson couplings to gauge bosons. Section III shows the
dominant channels for $pp\to HX\to VV^{*}X$ and $pp\to HHX\to
VV^{*}V^{\prime}V^{\prime*}X$ production events at the LHC and highlights the
sup-processes where the new heavy Higgs bosons can substantially modify these
cross-sections. This section also includes the calculation setup for our
analysis. In Section IV, we discuss our numerical results for single and
double Higgs production events in the MS-2HDM. Particularly, we show that the
cross-section of the $pp\to HHX\to VV^{*}V^{\prime}V^{\prime*}X$ processes is
significantly enhanced with respect to the SM. Finally, Section V contains our
conclusions.

##  II Type-II 2HDM and SM Alignment

The 2HDM contains two scalar iso-doublets, $\Phi_{i}\,(i=1,2)$, with
$\mathrm{U(1)_{Y}}$ hypercharges ${Y_{\Phi_{i}}=+1/2}$, as

${\Phi}_{i}\ =\ \left(\begin{matrix}\Phi_{i}^{+}\\\
\Phi_{i}^{0}\end{matrix}\right)$ (1)

The most general renormalisable 2HDM potential in terms of these doublets may
be conveniently written as

$\displaystyle V$ $\displaystyle=$
$\displaystyle\mu_{1}^{2}(\Phi_{1}^{\dagger}\Phi_{1})+\mu_{2}^{2}(\Phi_{2}^{\dagger}\Phi_{2})-\Big{[}m_{12}^{2}(\Phi_{1}^{\dagger}\Phi_{2})\>+\>{\rm
H.c.}\Big{]}$ (2) $\displaystyle+$
$\displaystyle\lambda_{1}(\Phi_{1}^{\dagger}\Phi_{1})^{2}+\lambda_{2}(\Phi_{2}^{\dagger}\Phi_{2})^{2}+\lambda_{3}(\Phi_{1}^{\dagger}\Phi_{1})(\Phi_{2}^{\dagger}\Phi_{2})+\lambda_{4}(\Phi_{1}^{\dagger}\Phi_{2})(\Phi_{2}^{\dagger}\Phi_{1})$
$\displaystyle+$ $\displaystyle\bigg{[}\,{1\over
2}\lambda_{5}(\Phi_{1}^{\dagger}\Phi_{2})^{2}+\lambda_{6}(\Phi_{1}^{\dagger}\Phi_{1})(\Phi_{1}^{\dagger}\Phi_{2})+\lambda_{7}(\Phi_{1}^{\dagger}\Phi_{2})(\Phi_{2}^{\dagger}\Phi_{2})\>+\>{\rm
H.c.}\bigg{]}\;,$

where the mass terms $\mu^{2}_{1,2}$ and quartic couplings $\lambda_{1,2,3,4}$
are real parameters. Instead, the remaining mass term $m^{2}_{12}$ and the
quartic couplings $\lambda_{5,6,7}$ are complex. Out of these 14 theoretical
parameters, only 11 are physical, since 3 parameters can be transformed away
using a SU(2) reparametrisation of the Higgs doublets. Here, we restrict our
attention to CP conservation and to CP-conserving vacua. In the Type-II 2HDM,
both scalar doublets receive nonzero vacuum expectation values (VEVs), as
$\langle\Phi_{1}\rangle=(0,\,v_{1}/\sqrt{2})^{\mathsf{T}}$ and
$\langle\Phi_{2}\rangle=(0,\,v_{2}/\sqrt{2})^{\mathsf{T}}$. Following the
standard linear expansion of the scalar doublets about their VEVs, these can
be re-expressed as

$\Phi_{i}\ =\ \left(\begin{matrix}\phi_{i}^{+}\\\
{1\over\sqrt{2}}(v_{i}+\phi_{i}+i\phi_{i}^{0})\end{matrix}\right)\;.$ (3)

Accordingly, the minimization conditions resulting from the 2HDM potential in
(2) take on the following forms:

$\displaystyle\mu_{1}^{2}\ $ $\displaystyle=$ $\displaystyle\
m_{12}^{2}{v_{2}\over v_{1}}-{1\over
2}v_{1}^{2}\big{(}2\lambda_{1}+3\lambda_{6}{{v_{2}\over
v_{1}}}+\lambda_{345}{({v_{2}\over v_{1}})^{2}}+\lambda_{7}{({v_{2}\over
v_{1}})^{3}}\big{)}\,,$ (4) $\displaystyle\mu_{2}^{2}\ $ $\displaystyle=$
$\displaystyle\ m_{12}^{2}{v_{1}\over v_{2}}-{1\over
2}v_{1}^{2}\big{(}2\lambda_{2}+3\lambda_{6}{v_{1}\over
v_{2}}+\lambda_{345}({v_{1}\over v_{2}})^{2}+\lambda_{7}({v_{1}\over
v_{2}})^{3}\big{)}\,,$ (5)

where ${\lambda}_{345}\equiv{\lambda}_{3}+{\lambda}_{4}+{\lambda}_{5}$,
${\lambda}_{5,6,7}\equiv{\rm Re}\,{\lambda}_{5,6,7}$ and
${m}_{12}^{2}\equiv{\rm Re}\,{m}_{12}^{2}$. Also, the combination of $v_{1}$
and $v_{2}$ form the VEV of the SM doublet,
$v^{2}\equiv{(v_{1}^{2}+v_{2}^{2})}$ and their ratio is read
$t_{\beta}\equiv\tan\beta=v_{2}/v_{1}$. Therefore, the mixing in the CP-even,
CP-odd and charged scalar sectors can be individually governed by the same
mixing angle $\beta$ as

$\begin{pmatrix}H_{\rm SM}\\\
H_{p}\end{pmatrix}=\mathcal{O}(\beta)\begin{pmatrix}\phi_{1}\\\
\phi_{2}\end{pmatrix},\quad\begin{pmatrix}G^{0}\\\
a\end{pmatrix}=\mathcal{O}(\beta)\begin{pmatrix}\phi_{1}^{0}\\\
\phi_{2}^{0}\end{pmatrix},\quad\begin{pmatrix}G^{\pm}\\\
h^{\pm}\end{pmatrix}=\mathcal{O}(\beta)\begin{pmatrix}\phi_{1}^{\pm}\\\
\phi_{2}^{\pm}\end{pmatrix},$ (6)

where the rotation matrix may be defined as below,

$\mathcal{O}(\beta)=\begin{pmatrix}c_{\beta}&s_{\beta}\\\
-s_{\beta}&c_{\beta}\end{pmatrix},$ (7)

with $c_{\beta}\equiv\cos\beta$ and $s_{\beta}\equiv\sin\beta$.

After the spontaneous symmetry breaking, the standard $W^{\pm}$ and $Z$ bosons
acquire their masses from the three would-be Goldstone bosons
$(G^{\pm},G^{0})$ Glashow:1961tr ; Kita:2011yv . Thereafter, the model remains
with five physical scalar states: two CP-even scalars ($H$ and $h$), one CP-
odd scalar $(a)$ and two charged bosons ($h^{\pm}$). The masses of the
$h^{\pm}=-s_{\beta}\phi_{1}^{\pm}+c_{\beta}\phi_{2}^{\pm}$ and
$a=-s_{\beta}\phi_{1}^{0}+c_{\beta}\phi_{2}^{0}$ scalars are given by

$\displaystyle M_{h^{\pm}}^{2}\ $ $\displaystyle=$ $\displaystyle\
{m_{12}^{2}\over s_{\beta}c_{\beta}}-{v^{2}\over
2}(\lambda_{4}+\lambda_{5})+{v^{2}\over
2s_{\beta}c_{\beta}}(\lambda_{6}c_{\beta}^{2}+\lambda_{7}s_{\beta}^{2}),$
$\displaystyle M_{a}^{2}\ $ $\displaystyle=$ $\displaystyle\
M_{h^{\pm}}^{2}+{v^{2}\over 2}(\lambda_{4}-\lambda_{5}),$ (8)

Moreover, the masses of the physical states $H$ and $h$ may be obtained by
diagonalising the CP-even $\phi_{1}$ and $\phi_{2}$ mass matrix $M^{2}_{S}$,

$M_{S}^{2}=\left(\begin{matrix}A\,&C\\\ C\,&B\end{matrix}\right),$ (9)

which may be explicitly written in the following form

$\displaystyle A=$ $\displaystyle
M_{a}^{2}s_{\beta}^{2}+v^{2}\left(2\lambda_{1}c_{\beta}^{2}+\lambda_{5}s_{\beta}^{2}+2\lambda_{6}s_{\beta}c_{\beta}\right),$
$\displaystyle B=$ $\displaystyle
M_{a}^{2}c_{\beta}^{2}+v^{2}\left(2\lambda_{2}s_{\beta}^{2}+\lambda_{5}c_{\beta}^{2}+2\lambda_{7}s_{\beta}c_{\beta}\right),$
$\displaystyle C=$ $\displaystyle-
M_{a}^{2}s_{\beta}c_{\beta}+v^{2}\left((\lambda_{3}+\lambda_{4})s_{\beta}c_{\beta}+\lambda_{6}c_{\beta}^{2}+\lambda_{7}s_{\beta}^{2}\right).$

The physical states $h$ and $H$ can be given in terms of the mixing angles
$\alpha$, as

$\begin{pmatrix}H\\\
h\end{pmatrix}=\mathcal{O}(\alpha)\begin{pmatrix}\phi_{1}\\\
\phi_{2}\end{pmatrix}=\mathcal{O}(\alpha)\mathcal{O}(\beta)^{-1}\begin{pmatrix}H_{\rm
SM}\\\ H_{p}\end{pmatrix}.$ (10)

Hence, the SM Higgs field may now be identified by the following linear fields
combination,

$H_{\text{SM}}\ =\ \phi_{1}\cos\beta+\phi_{2}\sin\beta\ =\
H\cos(\beta-\alpha)+h\sin(\beta-\alpha)\;.$ (11)

Henceforth, the SM-normalised couplings of the CP-even $H$ and $h$ scalars to
the EW gauge bosons are given by

$\displaystyle g_{hVV}=\sin(\beta-\alpha)\;,\qquad
g_{HVV}=\cos(\beta-\alpha).$ (12)

Thereby, there are two ways to realise exact SM alignment limit; (i) SM-like
$h$ scenario with $\sin(\beta-\alpha)=1$ and (ii) SM-like $H$ scenario with
$\cos(\beta-\alpha)=1$. Here, we consider the second realization with
$\beta=\alpha$ and assuming the CP-even scalar partner $h$ can be either
lighter or heavier than the observed scalar resonance at the LHC. As it can be
seen in (11), with this assumption, the SM-like Higgs boson becomes aligned
with one of the neutral eigenstates. So, the CP-even mass matrix $M_{S}^{2}$
takes on the form

$\displaystyle\widehat{M}_{S}^{2}$ $\displaystyle=$
$\displaystyle\mathcal{O}(\beta)M_{S}^{2}\mathcal{O}(\beta)^{-1}=\left(\begin{matrix}\widehat{A}&\widehat{C}\\\
\widehat{C}&\widehat{B}\end{matrix}\right),$ (13)

with

$\displaystyle\widehat{A}$ $\displaystyle=$ $\displaystyle
2v^{2}\left[c_{\beta}^{4}\lambda_{1}+s_{\beta}^{2}c_{\beta}^{2}\lambda_{345}+s_{\beta}^{4}\lambda_{2}+2s_{\beta}c_{\beta}\left(c_{\beta}^{2}\lambda_{6}+s_{\beta}^{2}\lambda_{7}\right)\right],$
$\displaystyle\widehat{B}$ $\displaystyle=$ $\displaystyle
M_{a}^{2}+\lambda_{5}v^{2}+2v^{2}\left[s_{\beta}^{2}c_{\beta}^{2}\left(\lambda_{1}+\lambda_{2}-\lambda_{345}\right)-s_{\beta}c_{\beta}\left(c_{\beta}^{2}-s_{\beta}^{2}\right)\left(\lambda_{6}-\lambda_{7}\right)\right],$
(14) $\displaystyle\widehat{C}$ $\displaystyle=$ $\displaystyle
v^{2}\left[s_{\beta}^{3}c_{\beta}\left(2\lambda_{2}-\lambda_{345}\right)-c_{\beta}^{3}s_{\beta}\left(2\lambda_{1}-\lambda_{345}\right)+c_{\beta}^{2}\left(1-4s_{\beta}^{2}\right)\lambda_{6}+s_{\beta}^{2}\left(4c_{\beta}^{2}-1\right)\lambda_{7}\right].$

From the above relations, we may observe that the SM alignment limit,
$\beta=\alpha$, can be understood if either $\widehat{C}\to 0$ or (ii)
$M_{h^{\pm}}\\!\sim\\!M_{a}\gg v$. However, these can be obtained naturally by
imposing the symmetries of model to constraint the certain parameters required
for the SM alignment limit. The 2HDM potential contains 13 accidental
symmetries, which have been fully classified in Pilaftsis:2011ed ;
Darvishi:2019dbh ; Birch-Sykes:2020btk . Of these, eight restrict the quartic
couplings such that the alignment condition $\widehat{C}\to 0$ is met.
However, only for three symmetries exact alignment can be achieved naturally
without imposing any constraints on the values of $\tan\beta$, nor on the
bilinear mass terms $\mu^{2}_{1,2}$ and $m^{2}_{12}$ Pilaftsis:2016erj ;
Darvishi:2020teg . In the simplest scenario, dubbed the MS-2HDM, the SM
alignment can be naturally realised as a consequence of an accidental SO(5)
symmetry in the Higgs sector Pilaftsis:2011ed ; Dev:2014yca ; Dev:2017org ;
Darvishi:2019ltl . In the MS-2HDM, the SO(5) symmetry puts severe restrictions
on the allowed form of the kinematic parameters of the 2HDM potential in (2),

$\displaystyle\mu_{1}^{2}=\mu_{2}^{2}\,,\quad
m_{12}^{2}=0\,,\quad\lambda_{2}=\lambda_{1}\,,\quad\
\lambda_{3}=2\lambda_{1}\,,\quad\lambda_{4}=\lambda_{5}=\lambda_{6}=\lambda_{7}=0\;.$
(15)

These parameters produce one massive CP-even scalar with
$M_{H}^{2}=2\lambda_{2}v^{2}$. The other four physical scalars $(h,a,h^{\pm})$
become massless pseudo-Goldstone bosons, that would participate in several
Higgs decay channels which are inconsistent with observation. To avoid these,
we consider that only the soft SO(5)-breaking parameter Re $m^{2}_{12}$ is
non-zero Dev:2014yca ; Dev:2017org ; Darvishi:2019ltl . With this minimal
addition to the MS-2HDM potential, the scalar-boson masses to a very good
approximation are given by

$M_{H}^{2}=2\lambda_{2}v^{2},\qquad M_{h}^{2}\ =\ M_{a}^{2}\ =\
M_{h^{\pm}}^{2}\ =\ {\,m_{12}^{2}\over s_{\beta}c_{\beta}}\;.$ (16)

In addition to softly broken the SO(5) symmetry by the bilinear scalar mass
term $m^{2}_{12}$, this symmetry can be broken explicitly by RG effects due to
non-zero values of gauge coupling $g^{\prime}$ and Yukawa couplings. The SO(5)
symmetry may be realised at large RG scales $\mu_{X}$ in two discrete values
$\mu^{(1)}_{X}\sim 10^{11}$ GeV and $\mu^{(2)}_{X}\sim 10^{18}$ GeV by
considering the threshold scales $\mu_{thr}=M_{h}^{\pm}=500$ GeV. Nonetheless,
assuming larger values for the threshold scale e.g. $\mu_{thr}=M_{h}^{\pm}=$1
TeV, 10 TeV and 100 TeV the RG scales $\mu_{X}$ can reach to $\sim\,10^{20}$
GeV Darvishi:2019ltl . In this context, when this symmetry breaks due to RG
effect, giving rise to calculable non-zero values for misalignment predictions
of all Higgs boson couplings to SM particles.

In the EW scale, the $H$-boson couplings in terms of the light-to-heavy
scalar-mixing parameter may be expressed by
$\theta_{\mathcal{S}}\equiv\widehat{C}/\widehat{B}$. So, the approximate
analytic expressions may be given by Darvishi:2019ltl

$\displaystyle g_{HVV}$ $\displaystyle\simeq 1-{\theta_{\mathcal{S}}^{2}\over
2},$ (17a) $\displaystyle g_{hVV}$
$\displaystyle\simeq-{\theta_{\mathcal{S}}}\ =\ {v^{2}s_{\beta}c_{\beta}\over
M_{a}^{2}+v^{2}\lambda_{5}}\,\Big{[}\,c_{\beta}^{2}\left(2\lambda_{1}-\lambda_{345}\right)-s_{\beta}^{2}\left(2\lambda_{2}-\lambda_{345}\right)\Big{]}.$
(17b)

Given the narrow experimental limits on the deviation of $g_{HVV}$ from 1, one
must have the parameter $\theta_{\mathcal{S}}\ll 1$. In the SM alignment
limit, we have $g_{Huu}\,\to\,1$ and $g_{Hdd}\to$ 1\. To this extent, an
approximate analytic expressions for the $h$\- and $H$-boson couplings to up-
and down-type quarks are given by

$\displaystyle g_{Huu}$ $\displaystyle\simeq
1+{t^{-1}_{\beta}}\,{\theta_{\mathcal{S}}},\qquad\qquad g_{Hdd}\simeq
1-{\theta_{\mathcal{S}}}\,{t_{\beta}},$ (18) $\displaystyle g_{huu}$
$\displaystyle\simeq-{\theta_{\mathcal{S}}}+{t^{-1}_{\beta}},\qquad\qquad
g_{hdd}\simeq-{\theta_{\mathcal{S}}}-{t_{\beta}}.$

Also, the trilinear Higgs couplings can be defined via the light-to-heavy
scalar-mixing parameter $\theta_{\mathcal{S}}$, in the following forms

$\displaystyle\kappa_{HHH}$ $\displaystyle\simeq\ {3\over
2v}\Big{[}\widehat{A}+({3\over 2v}\widehat{A}-{\,m_{12}^{2}\over
s_{\beta}c_{\beta}})\theta_{\mathcal{S}}^{2}\Big{]},$ (19)
$\displaystyle\kappa_{hHH}$ $\displaystyle\simeq\ {\theta_{\mathcal{S}}\over
2v}\Big{[}-2\widehat{A}-\widehat{B}+2{\,m_{12}^{2}\over
s_{\beta}c_{\beta}}\Big{]},$ (20)

where the $\kappa_{hHH}$ coupling vanishes in the exact alignment limit and
$\kappa_{HHH}$ converges to Its SM counterpart.

Couplings |  | ATLAS |  | CMS |  | $\tan\beta=2$ |  | $\tan\beta=5$ |  | $\tan\beta=50$
---|---|---|---|---|---|---|---|---|---|---
$|g_{HZZ}^{\text{low-scale}}|$ |  | [0.86, 1.00] |  | [0.90, 1.00] |  | 0.9999 |  | 0.9999 |  | 0.9999
$|g_{HZZ}^{\text{high-scale}}|$ |  |  |  |  |  | 0.9981 |  | 0.9998 |  | 0.9999
$|g_{Htt}^{\text{low-scale}}|$ |  | $1.31^{+0.35}_{-0.33}$ |  | $1.45^{+0.42}_{-0.32}$ |  | 1.0049 |  | 1.0014 |  | 1.0000
$|g_{Htt}^{\text{high-scale}}|$ |  |  |  |  |  | 1.0987 |  | 1.0179 |  | 1.0001
$|g_{Hbb}^{\text{low-scale}}|$ |  | $0.49^{+0.26}_{-0.19}$ |  | $0.57^{+0.16}_{-0.16}$ |  | 0.9803 |  | 0.9649 |  | 0.9590
$|g_{Hbb}^{\text{high-scale}}|$ |  |  |  |  |  | 0.8810 |  | 0.9264 |  | 0.9427
$|g_{HHH}^{\text{low-scale}}|$ |  | [-5.1, 12.0] |  | [-11.8, 18.8] |  | 0.9968 |  | 0.9968 |  | 0.9968
$|g_{HHH}^{\text{high-scale}}|$ |  |  |  |  |  | 0.9394 |  | 0.9936 |  | 0.9968
Table 1: Predicted values of the SM-like Higgs boson couplings to the $Z$
boson and to top- and bottom-quarks in the MS-2HDM for both scenarios with
low- and high-scale quartic coupling unification, assuming $M_{h^{\pm}}=500$
GeV. The corresponding central values for these couplings from ATLAS and CMS
are also given, including their uncertainties Khachatryan:2016vau . The ratio
of the Higgs boson self-coupling to its SM value, $g_{HHH}$, is constrained at
95$\%$ confidence level in observation (expectation) to
$-5.0<g_{HHH}<12.0\,(-5.8<g_{HHH}<12.0)$ and $-11.8<g_{HHH}<18.8$
$(-7.1<g_{HHH}<13.6)$ by ATLAS and CMS, respectively Aad:2019uzh ;
Sirunyan:2018ayu .

The MS-2HDM is a minimal and very predictive extension of the SM governed by
only three additional parameters: the unification scale $\mu_{X}$, the charged
Higgs mass $M_{h^{\pm}}$ and $\tan\beta$, allowing one to determine the entire
Higgs sector of the model. Previously, we presented our MS-2HDM benchmarks in
terms of these input parameters, for our misalignment predictions of the SM-
like Higgs-boson couplings to the $W^{\pm}$ and $Z$ bosons Darvishi:2019ltl .
These benchmarks are given in Table 1. Having considered the running mass
$M_{H}(M_{h^{\pm}}=500)$ to obtain $M_{H}(m_{t})\sim 125$ GeV, the mass of
heavy neutral Higgses are $M_{h}\sim 498$ GeV and $M_{a}\sim 495$ GeV.

Here, we intend to calculate the rate of $W^{\pm}/Z$-pair and -quadruplet
productions through single and double SM-like Higgs boson decay modes based on
these benchmarks.

##  III $W^{\pm}/Z$-Quadruplet Production and Even generation Framework

The $W^{\pm}/Z$-pair and -quadruplet production events are amongst the largest
Higgs-tagged signatures at the current LHC energies. However, the best
processes for probing the sings of an extended Higgs sector are the
$W^{\pm}/Z$\- quadruplet productions through double Higgs decays, i.e. $pp\to
HH\to VV^{*}V^{\prime}V^{\prime*}$. This is since, as we have shown in Section
II, the SM-like Higgs boson in the MS-2HDM couples to the EW gauge bosons with
coupling strength exactly as that of the SM Higgs boson, while the other
neutral heavy states do not couple to them at all Darvishi:2019ltl . On the
other hand, because of the mixing between the light and the heavy states, the
$W^{\pm}/Z$\- quadruplet production through double Higgs decays may be
enhanced. Henceforth, we calculate these processes within the MS-2HDM, which
can include Higgs trilinear decay modes $H/h\to HH$. The prediction for the
relative trilinear Higgs self-coupling to its SM value, $g_{HHH}$, for
different values of $\tan\beta$ which are allowed by the observed (expected)
data are summarized in Table 1. The $g_{HHH}$ is constrained at 95$\%$
confidence level in observation (expectation) to $-5.0<g_{HHH}<12.0$
$(-5.8<g_{HHH}<12.0)$ and $-11.8<g_{HHH}<18.8\,(-7.1<g_{HHH}<13.6)$ by ATLAS
and CMS, respectively Aad:2019uzh ; Sirunyan:2018ayu . Additionally, the non-
zero values for misalignment predictions giving rise to the sensitive BSM
decay mode $h\to HH$. The observed upper limit on the resonant production
cross-section times the branching fraction of $h\to HH$ with ATLAS detector
ranges between $40$ pb and 6.1 pb, while the expected limit ranges between
17.6 pb and 4.4 pb, for a hypothetical resonance with a mass in the range of
260-500 GeV Aaboud:2018ewm . Our prediction for the resonant production
$\sigma(gg\to h)\times Br(h\to HH)$ for $\tan\beta=2,\,5$ and $50$ are ranged
between 4.3 pb and 4.0 pb.

Note that there are several interesting channels that can be considered when
looking for the signal of the MS-2HDM model at the LHC, namely single and
double Higgs production events through $H\to b\bar{b}$, $H\to W^{+}W^{-}$,
$H\to ZZ$ and $H\to\tau\bar{\tau}$ decay channels. Indeed, the $H\to b\bar{b}$
decay has the largest branching ratio of $\sim\%57$, but they are generally
difficult to observe due to their large backgrounds. On the other hand, the
$H\to W^{+}W^{-}$ and $H\to ZZ$ decays have the combined branching ratio of
$\sim\%$25 and are very likely to produce clean experimental signatures via
their leptonic decay channels $W^{+}W^{-}\to
l^{+}\nu_{l}+l^{\prime-}\nu_{l^{\prime}}$ and $ZZ\to
l^{+}l^{-}+l^{\prime+}l^{\prime-}$.

Before delving into the details of this section, we check the $W^{\pm}/Z$-pair
production events through single Higgs decays to ensure the MS-2HDM
predictions are aligned with those of the SM. The exclusive production of the
EW gauge boson pairs and quadruplets via single and double Higgs decays can be
given as

$\displaystyle p(p_{1})+p(p_{2})$ $\displaystyle\to H(p_{H})+X$
$\displaystyle\to W^{+}(p_{W^{+}})+W^{-}(p_{W^{-}})+X\to
l_{1}^{+}\nu_{l_{1}}+l_{2}^{-}\nu_{l_{2}}+X,$ (21) $\displaystyle
p(p_{1})+p(p_{2})$ $\displaystyle\to H(p_{H})+X$ $\displaystyle\to
Z(p_{1,Z})+Z(p_{2,Z})+X\to l_{1}^{+}l_{1}^{-}+l_{2}^{+}l_{2}^{-}+X,$ (22)
$\displaystyle p(p_{1})+p(p_{2})$ $\displaystyle\to H(p_{1,H})+H(p_{2,H})+X$
$\displaystyle\to
W^{+}(p_{1,W^{+}})+W^{-}(p_{1,W^{-}})+W^{+}(p_{2,W^{+}})+W^{-}(p_{2,W^{-}})+X,$
$\displaystyle\to
l_{1}^{+}\nu_{l_{1}}+l_{2}^{-}\nu_{l_{2}}+l_{3}^{+}\nu_{l_{3}}+l_{4}^{-}\nu_{l_{4}}+X,$
(23) $\displaystyle p(p_{1})+p(p_{2})$ $\displaystyle\to
H(p_{1,H})+H(p_{2,H})+X$ $\displaystyle\to
Z(p_{1,Z})+Z(p_{2,Z})+Z(p_{3,Z})+Z(p_{4,Z})+X,$ $\displaystyle\to
l_{1}^{+}l_{1}^{-}+l_{2}^{+}l_{2}^{-}+l_{3}^{+}l_{3}^{-}+l_{4}^{+}l_{4}^{-}+X.$
(24)

In the above processes, the Higgs boson production mechanism is dominated by
the gluon fusion channels $gg\to H$ and $gg\to HH$, which account for
$\sim$%95 of the Higgs production rate at the LHC222In Masouminia:2018dwj , it
has been shown that one can simply enhance the LO differential cross-sections
of Higgs production with the use of a factorization-scale-dependant K-factor
and forget about the higher-order and virtual corrections to these processes,
e.g. via the $W/Z$ Higgs-strahlung sub-processes Baglio:2010ae . Nevertheless,
for the sake of completeness, we are considering the full range of real and
virtual QCD and QED contributions to (21), (22), (23) and (24), using the
Matchbox merging functionality within Herwig 7. Darvishi:2016fwo ;
Masouminia:2018dwj . Figure 1 displays all dominant sub-processes up to two
jets that are mediated through the exchange of a heavy virtual top/bottom
quark. The parts (a) and (b) of Figure 1 show the dominant single and double
Higgs production channels, respectively. The part (c) showcases the
contributions of the light and heavy Higgs bosons through the trilinear Higgs
vertices, where the signature of new physics would rise. Also, the leptonic
decay channels $W^{+}W^{-}\to l^{+}\nu_{l}+l^{\prime-}\nu_{l^{\prime}}$ and
$ZZ\to l^{+}l^{-}+l^{\prime+}l^{\prime-}$ have been considered to ensure a
clean observable signature and to prevent the reconstruction of $H\to VV^{*}$
resonances.

In our calculations for the production rates of the $pp\to HX\to VV^{*}$ and
$pp\to HHX\to VV^{*}V^{\prime}V^{\prime*}$ events, we utilize the Herwig 7
(v7.2.1) event generator Bahr:2008pv ; Bellm:2015jjp ; Bellm:2017bvx ;
Bellm:2019zci . This will be done for both the SM and the MS-2HDM, in the
collinear factorization framework. The contributing matrix elements are
generated with MadGraph5 and convoluted by MMHT2014 parton distribution
function libraries Harland-Lang:2014zoa  via LHAPDF interface Buckley:2014ana
. The underlying events are enhanced by an AO QCD+QED+EW parton shower scheme
Platzer:2009jq ; Platzer:2011bc ; Masouminia:2021  and the cluster model
hadronization Webber:1983if . The generated events are then analysed with
Rivet, based on the existing analysis MC_WWINC and MC_ZZINC which are modified
according to our needs.

Figure 1:  Dominant sub-processes for $W^{\pm}/Z$-pair and -quadruplet
productions through single and double Higgs boson decays at the LHC up to 2
jets; (a) $W^{\pm}/Z$-pair productions via single Higgs decay channels. (b)
$W^{\pm}/Z$-quadruplet productions through double Higgs boson decays. (c)
Double Higgs production channels via $H/h\to HH$ decays.

##  IV Numerical results and discussion

In this section, we present our results for $W^{\pm}/Z$-pair and quadruplet
production through the single and the double Higgs bosons decays at
$\sqrt{s}=13$ TeV. In the first step, we check the efficiency of our
calculation setup by evaluating the single Higgs production results in the SM
with the existing experimental observations from the CMS collaboration,
including their statistical and systematic uncertainties Darvishi:2019uzp ;
Sirunyan:2020tzo ; Sirunyan:2018sgc .

In Figure 2, we exhibit the results of our analysis for single Higgs bosons
production via $H\to VV^{*}$ decay channels. The top panel demonstrates the
kinematically reconstructed transverse momentum distribution of the exchanged
Higgs boson from the $H\to W^{+}W^{-}$ decay mode while the bottom panels
correspond to the transverse momentum and pseudo-rapidity from the $H\to ZZ$
decay mode. Figure 3 shows the fiducial cross-sections of single Higgs
production through $W^{\pm}$-pairs (left panel) and $Z$-pairs (right panel).
The event selection criteria for these calculations have been chosen in
accordance with the reported conditions in Sirunyan:2020tzo ; Sirunyan:2018sgc
. Despite the low precision of the experimental data, in both figures one can
readily observe that these predictions are perfectly capable of describing the
experimental measurements.

Figure 2: Differential cross-section for single Higgs boson production as a
function of its transverse momentum and pseudo-rapidity. The calculations are
in the colliner framewok, using Herwig 7 (v7.2.1) at $\sqrt{s}=13$ TeV. The
top panel corresponds to the $pp\to HX\to W^{+}W^{-}X$ channel while the
bottom panels are for the $pp\to HX\to ZZX$ channel. The data are from the CMS
collaboration Sirunyan:2020tzo ; Sirunyan:2018sgc . To calculate the
uncertainty region, we have manipulated the factorization hard-scale by a
factor of 2. Figure 3: Feducial cross-section for $W^{\pm}(Z)$-pair
productions through single SM Higgs boson decays at the LHC for $\sqrt{s}=13$
TeV displayed in left (right) panel. The data are from the CMS collaboration,
including their statistical and systematic uncertainties Sirunyan:2020tzo ;
Sirunyan:2018sgc . To calculate the uncertainty region, we have manipulated
the factorization hard-scale by a factor of 2.

In the next step, we calculate the single and the double Higgs production
events within the MS-2HDM. According to our discussion in Section II, the
MS-2HDM has two conformally-invariant quartic coupling unification points
$\mu^{(1)}_{X}$ (low-scale) and $\mu^{(2)}_{X}$ (high-scale), for a given
choice of the charged Higgs-boson mass $M_{h^{\pm}}$ and $\tan\beta$
Darvishi:2019ltl ; Darvishi:2020teg . Thus, we perform these analysis for both
unification points. In Figures 4 and 5, we display the differential cross-
sections for $pp\to H\to VV^{*}$ production as a function of pseudo-rapidity
of the produced gauge bosons in the MS-2HDM. The results for the low-scale
(LS) and the high-scale (HS) points are shown in the left and the right
panels, respectively. The top panels correspond to the kinematic properties of
the EW bosons while the bottom panels depict the behaviour of gauge boson
pairs. The calculations have been done at $13$ TeV center-of-mass energy for
$M_{h^{\pm}}=500~{}\rm GeV$ and the typical values of $\tan\beta$, such as
$\tan\beta=2,\,5$ and $50$, relevant to the benchmarks of Table 1. By analogy,
Figures 6 and 7 show the transverse momentum distributions for $W^{\pm}$ and
$Z$ bosons production through single SM-like Higgs bosons.

As expected, we observe that the results of $pp\to H\to W^{+}W^{-}$ and $pp\to
H\to ZZ$ for both lower- and higher-scale quartic coupling unification points
are in excellent agreement with the SM and the experimental data. Obviously,
this is since the normalised couplings $g_{HVV}$ approaches the SM value
$g_{H_{\text{SM}}VV}=1$ for both points. However, the deviation of the MS-2HDM
results from the SM predictions is higher because the degree of misalignment
reaches its maximum value for $\tan\beta=2$, while still remaining within
their $1\sigma$ uncertainty.

Figure 4: Differential cross-section as a function of pseudo-rapidity for
$W^{\pm}$-pair productions through single SM-like Higgs boson ($H$) events
displayed for lower-scale (higher-scale) quartic coupling unification point in
left panels (right panels). These are shown for different values of
$\tan\beta$ at $\sqrt{s}=13$ TeV. The top panels correspond to the kinematic
properties of the $W^{\pm}$ bosons while the bottom panels depict the
behaviour of $W^{\pm}$-pairs. Note that, these results are compared with
relevant theoretical predictions in the SM within its uncertainty bounds.

Figure 5: The same as in Figure 4, but for $Z$-pair productions

.

Figure 6: The same as in Figure 4, but for differential cross-section as a
function of transverse momentum for $W^{\pm}$-pair productions. Figure 7: The
same as in Figure 4, but for differential cross-section as a function of
transverse momentum for $Z$-pair productions.

Now, let us turn our attention to the $W^{\pm}/Z$-quadruplet production events
through double Higgs decay channels, i.e. $pp\to HH\to
VV^{*}V^{\prime}V^{\prime*}$. Here, the sub-processes involving $H/h\to HH$
decays may have large contributions into the $pp\to HH$ production rate, as
shown in Figure 1(c). In Figures 8 and 9, we exhibit our results for the
production of these events with $\tan\beta=2$, $5$ and $50$ at $13$ TeV
center-of-mass energy. In both figures, the top panels represent the
reconstructed kinematics of the exchanged $H$ scalars while bottom plots show
the rates of production as functions of the pseudo-rapidity of the gauge
bosons. These are shown for both the low-scale and the high-scale unification
points. From these plots, we observe a substantial increase in the production
rate of the SM-like Higgs bosons compared to their SM counterparts. This
becomes more pronounced for the smaller values of $\tan\beta$, which has a
nearly 2-fold increase compared to similar SM prediction. Despite the fact
that these processes have expectedly smaller cross-sections in comparison with
the single Higgs production events, they have substantial deviance from the SM
and may be directly observed at the LHC data.

In a similar fashion, Figures 10 and 11 demonstrate the transverse momentum
distributions of the differential cross-section for the double Higgs
production events. Observe that, the peaks in the Higgs and $W^{\pm}/Z$ bosons
$p_{\perp}$ distributions are increased by a factor $\sim$3 with respect to
the SM for $p_{\perp}<200$ GeV. However, the $p_{\perp}$ distribution’s tails
converge to the SM predictions in the high-$p_{\perp}$ regions. Therefore, the
signature of the MS-2HDM may be observed in the low-$p_{\perp}$ regions of the
$W^{\pm}/Z$-quadruplet production events through double Higgs decay channels.
Our observations can be readily generalised to other realisations of the 2HDM
in their alignment limits.

Figure 8: Differential cross-section as a function of pseudo-rapidity for
$W^{\pm}$-quadruplet productions through double SM-like Higgs boson ($H$)
events displayed for lower-scale (higher-scale) quartic coupling unification
point in left panels (right panels). These are shown for different values of
$\tan\beta$ at $\sqrt{s}=13$ TeV. The top panels correspond to the kinematic
properties of the $W^{\pm}$ bosons while the bottom panels depict the
behaviour of $W^{\pm}$-pairs. Note that, these results are compared with
relevant theoretical predictions in the SM within its uncertainty bounds.
Figure 9: The same as in Figure 8, but for $Z$-quadruplet productions. Figure
10: The same as in Figure 8, but for differential cross-section as a function
of transverse momentum for $W^{\pm}$-quadruplet productions. Figure 11: The
same as in Figure 8, but for differential cross-section as a function of
transverse momentum for $Z$-quadruplet productions.

As a final remark, it is also possible to look for the signature of the
MS-2HDM through other significant signals like $b\bar{b}b\bar{b}$ and
$W^{+}W^{-}b\bar{b}$ decay channel. In Figure 12, we show the differential
cross-section for the production of Higgs pairs with the subsequent $HH\to
W^{+}W^{-}b\bar{b}$ decays. It is noticeable that the resulting signals are
expectedly larger compared to their counterparts in the $HH\to
VV^{*}V^{\prime}V^{\prime*}$ channels. This can be contributed to the larger
branching ratio of $Br(H\to b\bar{b})\sim\%57$ compared to $Br(H\to
W^{+}W^{-}/ZZ)\sim\%$25 as well as direct contributions of the heavy neutral
Higgs bosons through $h\to b\bar{b}$ ($a\to b\bar{b}$) decay mode with
$Br(h\to b\bar{b})\sim\%$5 ($Br(a\to b\bar{b})\sim\%$2). However, these
signals despite producing numerically larger contributions are plagued with
very large SM backgrounds. In contrast, the $VV^{*}V^{\prime}V^{\prime*}$
decay channels with their leptonic final states produce large and clean
signatures for the MS-2HDM that might be observed in the out-coming LHC data.

Figure 12: Differential cross-section for the production of Higgs pairs with
the subsequent $HH\to W^{+}W^{-}b\bar{b}$ decays. The differential cross-
sections are plotted as functions of pseudo-rapidity (top panels) and
transverse momentum (bottom panels) of the $W^{\pm}$ bosons (left panels) and
$b$-quarks (right panels) productions for higher-scale quartic coupling
unification point. These are shown for different values of $\tan\beta$ at
$\sqrt{s}=13$ TeV. Note that, these results are compared with relevant
theoretical predictions in the SM within its uncertainty bounds.

##  V Conclusions

The data collected from CERN’s LHC impose constraints over the coupling
strengths of the Higgs boson, primarily to the EW gauge bosons
($V=W^{\pm},Z$), which are very close to the SM predictions. This simple fact
severely restricts the form of possible scalar-sector extensions of the SM. In
this study, we have considered the MS-2HDM where the SM alignment can be
achieved naturally by the virtue of an SO(5) symmetry imposed on the 2HDM. The
MS-2HDM is a minimal and very predictive extension of the SM governed by only
three parameters: the unification scale $\mu_{X}$, the charged Higgs mass
$M_{h^{\pm}}$ and $\tan\beta$ which allow one to determine the entire Higgs
sector of the model.

Given the remarkable features of the MS-2HDM Darvishi:2019ltl ;
Darvishi:2020teg , we have investigated the possible signature of this model
via $W^{\pm}/Z$-quadruplet productions at the LHC. We have performed our
calculations with NLO QCD accuracy for $pp\to HX\to VV^{*}X$ and $pp\to HHX\to
VV^{*}V^{\prime}V^{\prime*}X$ processes for different values of $\tan\beta$,
using the Herwig 7 multi-purpose event generator at $\sqrt{s}=13$ TeV center-
of-mass energy. The corresponding amplitudes are provided by MadGraph5, up to
one QCD loop and two jets. The produced underlying events are showered by an
AO MC@NLO matched QCD+QED+EW parton shower and the results have been analysed
using Rivet.

We have shown that the predictions for $W^{\pm}/Z$-pair productions through
single SM-like Higgs boson events are aligned with their SM counterparts,
while the presence of the heavy Higgs states significantly enhances the
$W^{\pm}/Z$-quadruple productions. Particularly, we have found that the cross-
section for these events is increased by a factor $\sim$3 in $p_{\perp}<200$
GeV region with respect to the SM. However, these distributions converge to
the SM predictions in the high-$p_{\perp}$ regions. Therefore, the signature
of the MS-2HDM may be observed in the low-$p_{\perp}$ regions of the
$W^{\pm}/Z$-quadruplet production events through double Higgs decay channels.
These observations are very helpful toward the possible future discovery of
this model.

###### Acknowledgements.

The authors would like to thank Prof. A. Pilaftsis and Prof. P. Richardson for
their instructive discussions. ND is supported by the Lancaster-Manchester-
Sheffield Consortium for Fundamental Physics, under STFC research grant
ST/P000800/1. The work of ND is also supported in part by the National Science
Centre (Poland) as a research project, decision no 2017/25/B/ST2/00191 and by
the Polish National Science Centre HARMONIA grant under contract
UMO-2015/20/M/ST2/00518 (2016-2020). MRM is supported by the UK Science and
Technology Facilities Council (grant numbers ST/P001246/1). This work has
received funding from the European Union’s Horizon 2020 research and
innovation program as part of the Marie Skłodowska-Curie Innovative Training
Network MCnetITN3 (grant agreement no. 722104).

## References

  * (1) S. L. Glashow, Nucl. Phys. 22 (1961), 579-588
  * (2) J. Goldstone, A. Salam and S. Weinberg, Phys. Rev. 127 (1962), 965-970
  * (3) S. Weinberg, Phys. Rev. Lett. 19 (1967), 1264-1266
  * (4) A. Salam, Elementary Particle Physics: Relativistic Groups and Analyticity, Eighth Nobel Symposium (Almqvist and Wiksell, Stockholm, 1968), p. 367.
  * (5) S. L. Glashow, Nucl. Phys. 22 (1961), 579-588.
  * (6) F. Englert and R. Brout, Phys. Rev. Lett. 13 (1964), 321-323
  * (7) P. W. Higgs, Phys. Rev. Lett. 13 (1964), 508-509
  * (8) [ATLAS], ATLAS-CONF-2019-004.
  * (9) [CMS], CMS-PAS-SMP-19-001.
  * (10) I. F. Ginzburg, M. Krawczyk and P. Osland, [arXiv:hep-ph/9909455 [hep-ph]].
  * (11) A. Delgado, G. Nardini and M. Quiros, JHEP 07 (2013), 054.
  * (12) P. H. Chankowski, T. Farris, B. Grzadkowski, J. F. Gunion, J. Kalinowski and M. Krawczyk, [arXiv:hep-ph/0009271 [hep-ph]].
  * (13) M. Carena, I. Low, N. R. Shah and C. E. M. Wagner, JHEP 04 (2014), 015.
  * (14) P. S. Bhupal Dev and A. Pilaftsis, JHEP 12 (2014), 024 [erratum: JHEP 11 (2015), 147].
  * (15) J. Bernon, J. F. Gunion, H. E. Haber, Y. Jiang and S. Kraml, Phys. Rev. D 92 (2015) no.7, 075004.
  * (16) K. Benakli, M. D. Goodsell and S. L. Williamson, Eur. Phys. J. C 78 (2018) no.8, 658.
  * (17) K. Lane and W. Shepherd, Phys. Rev. D 99 (2019) no.5, 055015.
  * (18) N. Darvishi, JHEP 11 (2016), 065.
  * (19) N. Darvishi and M. Krawczyk, Nucl. Phys. B 926 (2018), 167-178.
  * (20) T. D. Lee, Phys. Rev. D 8 (1973), 1226-1239
  * (21) A. Pilaftsis and C. E. M. Wagner, Nucl. Phys. B 553 (1999), 3-42
  * (22) G. C. Branco, M. N. Rebelo and J. I. Silva-Marcos, Phys. Lett. B 614 (2005), 187-194.
  * (23) V. Silveira and A. Zee, Phys. Lett. B 161 (1985), 136-140.
  * (24) C. Bonilla, D. Sokolowska, N. Darvishi, J. L. Diaz-Cruz and M. Krawczyk, J. Phys. G 43 (2016) no.6, 065001.
  * (25) M. Krawczyk, N. Darvishi and D. Sokolowska, Acta Phys. Polon. B 47 (2016), 183.
  * (26) V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett. B 155 (1985), 36.
  * (27) A. G. Cohen, D. B. Kaplan and A. E. Nelson, Ann. Rev. Nucl. Part. Sci. 43 (1993), 27-70.
  * (28) R. A. Battye, G. D. Brawn and A. Pilaftsis, JHEP 08 (2011), 020.
  * (29) A. Pilaftsis, Phys. Lett. B 706 (2012), 465-469.
  * (30) N. Darvishi and A. Pilaftsis, Phys. Rev. D 101 (2020) no.9, 095008.
  * (31) N. Darvishi and A. Pilaftsis, PoS CORFU2019 (2020), 064.
  * (32) N. Darvishi and A. Pilaftsis, Phys. Rev. D 99 (2019) no.11, 115014.
  * (33) A. Pilaftsis, Phys. Rev. D 93 (2016) no.7, 075012.
  * (34) P. S. Bhupal Dev and A. Pilaftsis, J. Phys. Conf. Ser. 873 (2017) no.1, 012008.
  * (35) E. Hanson, W. Klemm, R. Naranjo, Y. Peters and A. Pilaftsis, Phys. Rev. D 100 (2019) no.3, 035026.
  * (36) C. Birch-Sykes, N. Darvishi, Y. Peters and A. Pilaftsis, Nucl. Phys. B 960 (2020), 115171.
  * (37) N. Darvishi and M. R. Masouminia, Nucl. Phys. B 923 (2017), 491-507.
  * (38) N. Darvishi, M. R. Masouminia and K. Ostrolenk, Phys. Rev. D 101 (2020) no.1, 014007.
  * (39) G. Aad et al. [ATLAS], Phys. Lett. B 800 (2020), 135103 [arXiv:1906.02025 [hep-ex]].
  * (40) A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 122 (2019) no.12, 121803 [arXiv:1811.09689 [hep-ex]].
  * (41) M. Aaboud et al. [ATLAS], Eur. Phys. J. C 78 (2018) no.12, 1007 [arXiv:1807.08567 [hep-ex]].
  * (42) M. Bahr, S. Gieseke, M. A. Gigg, D. Grellscheid, K. Hamilton, O. Latunde-Dada, S. Platzer, P. Richardson, M. H. Seymour and A. Sherstnev, et al. Eur. Phys. J. C 58 (2008), 639-707.
  * (43) J. Bellm et al., Eur. Phys. J. C 76 (2016) no.4, 196.
  * (44) J. Bellm et al., arXiv:1705.06919 [hep-ph].
  * (45) J. Bellm, G. Bewick, S. Ferrario Ravasio, S. Gieseke, D. Grellscheid, P. Kirchgaeßer, M. R. Masouminia, G. Nail, A. Papaefstathiou and S. Platzer, et al. Eur. Phys. J. C 80 (2020) no.5, 452.
  * (46) J. Alwall et al., JHEP 1407 (2014) 079.
  * (47) J. Bellm, C. B. Duncan, S. Gieseke, M. Myska and A. Siódmok, Eur. Phys. J. C 79 (2019) no.12, 1003.
  * (48) S. Platzer and S. Gieseke, Eur. Phys. J. C 72 (2012), 2187.
  * (49) S. Frixione and B. R. Webber, JHEP 0206 (2002) 029.
  * (50) M. R. Masouminia, P. Richardson, [to appear soon].
  * (51) A. Buckley, J. Butterworth, L. Lonnblad, D. Grellscheid, H. Hoeth, J. Monk, H. Schulz and F. Siegert, Comput. Phys. Commun. 184 (2013), 2803-2819.
  * (52) T. Kita, J. Phys. Soc. Jap. 80 (2011), 084606.
  * (53) G. Aad et al. [ATLAS and CMS], JHEP 08 (2016), 045.
  * (54) M. Modarres, M. R. Masouminia, R. Aminzadeh Nik, H. Hosseinkhani and N. Olanj, Nucl. Phys. B 926 (2018), 406-426.
  * (55) J. Baglio and A. Djouadi, JHEP 03 (2011), 055.
  * (56) L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, Eur.Phys.J.C 75(2015)204.
  * (57) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr and G. Watt, Eur. Phys. J. C 75 (2015), 132.
  * (58) S. Platzer and S. Gieseke, JHEP 01 (2011), 024.
  * (59) B. R. Webber, Nucl. Phys. B 238 (1984), 492-528.
  * (60) A. M. Sirunyan et al. [CMS], [arXiv:2007.01984 [hep-ex]].
  * (61) A. M. Sirunyan et al. [CMS], Phys. Lett. B 792 (2019), 369-396.


# HyperGALE: ASD Classification via Hypergraph Gated Attention with Learnable
HyperEdges

Mehul Arora, Chirag S. Jain, Lalith Bharadwaj Baru, Kamalaker Dadi, Bapi Raju
S  Brain Cognition Computation Lab
IIIT Hyderabad, India
mehul.arora<EMAIL_ADDRESS>
###### Abstract

Autism Spectrum Disorder (ASD) is a neurodevelopmental condition characterized
by varied social cognitive challenges and repetitive behavioral patterns.
Identifying reliable brain imaging-based biomarkers for ASD has been a
persistent challenge due to the spectrum’s diverse symptomatology. Existing
baselines in the field have made significant strides in this direction, yet
there remains room for improvement in both performance and interpretability.
We propose _HyperGALE_ , which builds upon the hypergraph by incorporating
learned hyperedges and gated attention mechanisms. This approach has led to
substantial improvements in the model’s ability to interpret complex brain
graph data, offering deeper insights into ASD biomarker characterization.
Evaluated on the extensive ABIDE II dataset, _HyperGALE_ not only improves
interpretability but also demonstrates statistically significant enhancements
in key performance metrics compared to both previous baselines and the
foundational hypergraph model. The advancement _HyperGALE_ brings to ASD
research highlights the potential of sophisticated graph-based techniques in
neurodevelopmental studies. The source code and implementation instructions
are available at Github.

###### Index Terms:

ASD Classification, Hypergraphs, Graph Neural Networks, ABIDE II.

##  I Introduction

Figure 1: ASD Processing Pipeline and the Proposed HyperGale Architecuture.
(a) fMRI data is preprocessed with Schaefer’s parcellation and later converted
into a functional connectivity (FC) matrix. (b) The FC matrix is then
converted into a hypergraph with learnable hyperedge weights. Subsequently,
the hypergraph features are sent to a series of hypergraph convolution layers
and gated attention is applied in order to capture the importance of relevant
information. Finally, the information extracted after attention layers is
aggregated using the readout function. Finally, from these lower dimensional
features of readout, a sigmoid activation leads to binary classification of
ASD vs Typically Developing (TD).

Autism Spectrum Disorder (ASD) is a complex neurodevelopmental condition
characterized by challenges in social cognition and repetitive behaviors.
Affecting a significant portion of the population, with a prevalence of
approximately 1 in 36 children [1], ASD’s diagnosis is complicated by its
heterogeneous nature and the limitations of existing diagnostic criteria.
Emerging research suggests that the atypical behaviors observed in ASD may be
linked to distinct patterns in functional connectivity (FC) within the brain
[2][3].

Functional connectivity represents dynamic interactions and synchronizations
between various brain regions. These interactions form intricate networks that
are foundational for cognitive processes and behaviors, exemplified by
networks such as the Default Mode Network (DMN) and Dorsal Attention Network
(DAN) [4]. The importance of these networks in ASD is increasingly recognized
in contemporary research [5], underscoring the need for analytical models
capable of discerning complex FC patterns associated with ASD.

Several models have been applied to understand ASD, ranging from traditional
methods like SVM and Random Forest [6] to neural networks, CNNs and
Transformers [7, 8, 9, 10, 11]. Graph-based methods have also been explored,
retaining complex brain information for ASD classification. However, these
methods often overlook higher-order relationships in the brain’s network,
focusing on dyadic interactions. This limitation is significant in the context
of ASD, where the implicated brain regions can be highly variable and diverse
[12].

Addressing the intricate challenges in ASD diagnosis, our research introduces
_HyperGALE_ , at the confluence of computational neuroscience and graph
theory. Unlike traditional graph-based methods that primarily focus on
pairwise interactions, _HyperGALE_ utilizes hypergraph convolutions. This
approach allows the model to capture high-order relationships within the
brain’s network, crucial for understanding the complex and heterogeneous
nature of ASD. Coupled with gated attention, _HyperGALE_ not only discerns
intricate patterns in functional connectivity but also provides a complex
understanding of brain region interactions associated with ASD. Our key
contributions are as follows:

  1. 1.

We developed _HyperGALE_ , a novel ASD classification approach utilizing
modified hypergraph convolution and gated attention to identify critical brain
regions implicated in ASD. Additionally, we re-implemented and compared
several holistic baseline models on the ABIDE-II dataset. These comparisons
are based on accuracy, AUC, sensitivity, and specificity, highlighting
_HyperGALE_ ’s superior performance.

  2. 2.

Our investigation into the model’s hyperparameters, such as the number of
Regions of Interest (ROIs) in hyperedges and the hypergraph layer count,
provided insights into their impact on _HyperGALE_ ’s performance. We also
assessed the robustness of all models, including _HyperGALE_ , through
multiple runs with different initializations and dataset distributions, noting
_HyperGALE_ ’s consistently low standard deviation in outcomes. The model’s
generalization capabilities were further evidenced by employing a leave-one-
site-out strategy.

  3. 3.

The advancements in _HyperGALE_ go beyond performance metrics, offering
interpretative insights into ASD’s qualitative aspects through learnt
hyperedges and gated attention. This bridges the gap between computational
analysis and neuroscience, providing a deeper understanding of
neurodevelopmental disorders such as ASD.

##  II Related Works

In this section, we begin by briefly reviewing traditional machine-learning,
deep-learning architectures that do not incorporate graph structures. We then
transition to discussing the evolution and significance of graph-based and
hypergraph-based methods in the context of ASD classification.

Figure 2:  Modified hypergraph convolution proposed in HyperGALE pipeline.
Starting with ROI features derived from the fully connected matrix, these
initial inputs undergo node-to-hyperedge propagation to form hyperedge
features (step (b)). These features are scaled using learned weights, as in
step (c), the activation happens in different hyperedges with different
values. Finally, the scaled hyperedge features are propagated to new node
features (step (d)).

###  II-A Overview of Non-Graph Methods

Traditional machine learning methods, such as Support Vector Machines (SVM),
Random Forests, and Gradient Boosting Classifiers (GBC)[6], have been applied
to ASD classification. These methods, however, often face challenges in
handling the high variance inherent in ASD data, attributed to site-specific
variations and the heterogeneous nature of the disorder.

In the realm of deep learning, architectures like CNNs and Transformers have
shown promise. Eslami _et al._ [8] implemented an autoencoder with data
augmentation, while Kawahara _et al._ [9] utilized a CNN, treating the
adjacency matrix as an image rather than exploiting its inherent graph
structure. Kan _et al._ [11] innovatively employed a transformer-based
approach for node embeddings with self-supervised clustering for readout.

###  II-B Graph-based Methods

A myriad of graph-based models was proposed on top of the standard GNNs [13],
[14], [15]. For example, Cao _et al._ [16] constructed a deep ASD diagnostic
framework based on 16-layer GCN with ResNet units and DropEdge strategy to
avoid certain problems such as vanishing gradient, and over-smoothing. Kazi
_et al._ [17] proposes ‘inception modules’ which are capable of capturing
intra- and inter-graph structural heterogeneity during convolutions. Yao _et
al._ [18] introduce a multi-scale triplet graph convolutional network that
employs multi-scale templates and a triplet GCN (TGCN) model to learn multi-
scale graph representations of brain FC. A recent significant study by Chen
_et al._ [7] put forth a competitive model that used both the resting state
fMRI (brain function) and T1 weighted MRI (brain structure) to generate
separate node and edge embeddings using a Transformer block.

Since the FC matrix is symmetric, traditional models tend to use the
upper/lower triangular matrix features of FC, thereby lacking message passing.
Therefore, these methods struggle to aggregate local and global information,
resulting in poor performance. The CNNs benefit from translational
equivariance and the Transformer’s key feature is the self-attention
mechanism. Although, these methods have good properties, due to a lack of
message passing they fail to capture aggregated neighbourhood information.
Graph-based approaches on FC matrices are applied to understand interactions
among regions of interest, but current methods can not take advantage of
higher-order proximity information due to reliance on static dyadic edges. To
address these limitations and accommodate the heterogeneity in individuals
with ASD, exploring hypergraph models is considered to be a promising step.

###  II-C Hypergraph-based Methods

Recent ASD classification studies have explored hypergraph frameworks, such as
Hypergraph U-Net[19] and multi-view HGNN[20], utilizing limited multi-modal
data. Frameworks like those of Shao _et al._ [21] for Alzheimer’s
classification incorporate feature selection and group-sparsity regularization
in their hypergraph approach. Similarly, Liu _et al._ [22] use a view-aligned
regularizer for multimodal coherence in hypergraph creation. Another notable
work, [23], employs dynamic hypergraph neural networks, utilizing KNN and
KMeans for iterative hypergraph generation in ASD classification.

Different from these methods, _HyperGALE_ utilizes the brain’s inherent
hypergraph structure with a specialized hypergraph convolution technique and
learnable hyperedges. This approach enhances the interpretation of complex
brain networks for more effective ASD classification.

##  III Methodology

In this section, we initially discuss each component of our method step-by-
step and eventually detail the entire pipeline.

###  III-A Dataset and Preprocessing

To systematically benchmark various graph learning techniques for fMRI-based
disease classification, we used large-scale openly-accessible fMRI dataset
from ABIDE-II
consortium111http://fcon_1000.projects.nitrc.org/indi/abide/abide_II.html.
This work included 812 subjects, out of which ASD: 384 and TD: 428,
respectively, accumulated across 16 sites as shown on the left of Figure 4.
The resting-state fMRI dataset was pre-processed using the default processing
pipeline considering global signal regression and bandpass filtering
implemented with C-PAC pipeline [24].

Briefly, the pre-processing pipeline included: removing the skull regions and
segmentation of each anatomical image into three tissues, followed by
normalization to the MNI152 template using Advanced Normalization Tools (ANTs)
[25]. fMRI pre-processing included slice timing correction, motion correction,
global mean intensity normalization as well as nuisance signal regression. In
nuisance signal regression, the number of parameters used were 24 parameters
for head motion, CompCor [26] with five principal components signals from
Cerebrospinal fluid and white matter, linear and quadratic trends for low-
frequency drifts. fMRI images were then co-registered with subject-specific
corresponding anatomical images and finally normalized to the MNI152 space
using ANTs.

After pre-processing, timeseries signals were extracted for each node where
each node belongs to 7 networks 400 parcellations from Schaefer atlas [27]. On
each node time series, functional connectivity estimation using the Ledoit-
Wolf regularized shrinkage estimator [28] and full correlation was implemented
with Nilearn package [29].

###  III-B Hypergraph Modeling

As introduced earlier, we construct a graph based on the FC correlation
matrix. Also, the fundamental notations and the formulations are adopted from
Feng _et al._ [30]. Basic Notations: Suppose
$\mathcal{G}=(\mathcal{V},\mathcal{E})$ be our graph with $\mathcal{V}$ as
vertices and $\mathcal{E}$ edges respectively. Where,
$\mathcal{V}\in\\{v_{1},v_{2},\cdots v_{N}\\}$ and
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$. The Adjacency matrix
$\mathbb{A}\in\mathbb{R}^{N\times N}$ informs us of the pairwise interactions
between each and every node.

Similarly, to construct a hypergraph
$\mathcal{G}_{hyper}=(\mathcal{V},\mathcal{E}_{h})$ with $N$ nodes and $K$
hyperedges. Similar to the adjacency matrix in the graph, in hypergraph we
have _Incidence matrix_ $\mathcal{H}\in\mathbb{R}^{N\times K}$. Here, the
hyperedges $e_{j}\in\mathcal{E}_{h}$ ($j=1,2,\cdots K$) are dumped in a
diagonal weight matrix $W\in\mathbb{R}^{K\times K}$. The entries in the
aforementioned incidence matrix can be represented as

$\mathcal{H}_{ij}=\begin{cases}1&\quad\text{if }v_{i}\in e_{j}\\\
0&\quad\text{Else.}\\\ \end{cases}$

The vertex degree and the hyperedge degree of the hypergraph are defined as
$\mathcal{D}_{ii}=\sum_{j=1}^{K}W_{jj}\mathcal{H}_{ij}$ and
$\mathcal{B}_{jj}=\sum_{i=1}^{N}\mathcal{H}_{ij}$ respectively.

TABLE I: The Table illustrates the performance of our proposed method
HyperGALE as compared to the existing state-of-the-art approaches. The train
and test proportions considered for the study are 90% and 10%, respectively.
The experiments are conducted for 5 different non-overlapping proportions of
data and the mean and standard deviations are reported.

Modality | Methods | Performance Metrics
---|---|---
Accuracy | AUC | Sensitivity | Specificity
Traditional Methods | SVM [31] | 68.83$\pm$ 2.44 | 68.07$\pm$ 2.49 | 78.82$\pm$ 2.65 | 57.31$\pm$ 3.01
Random Forest [32] | 60.00$\pm$ 5.55 | 58.80$\pm$6.11 | 70.81$\pm$ 4.11 | 48.90$\pm$7.24
Gradient Boosting | 62.01$\pm$ 2.99 | 61.45$\pm$ 2.92 | 68.86$\pm$ 4.54 | 54.05$\pm$ 5.23
Non-graph Methods | ASD-Diagnet [8] | 70.8 $\pm$1.49 | 71.12 $\pm$2.06 | 68.36 $\pm$8.13 | 72.2 $\pm$10.56
BrainNetCNN [9] | 67.32$\pm$ 4.24 | 69.07$\pm$ 3.05 | 68.89$\pm$ 8.62 | 65.40$\pm$ 11.56
Graph based Methods | GCN [13] | 72.68$\pm$3.49 | 75.65$\pm$ 3.57 | 79.56$\pm$ 8.42 | 64.33$\pm$ 9.65
GAT [14] | 68.53$\pm$ 3.40 | 71.70$\pm$ 4.11 | 65.78$\pm$ 9.76 | 71.89$\pm$ 9.75
GraphSAGE [15] | 71.22$\pm$2.55 | 74.64$\pm$4.60 | 75.43$\pm$17.88 | 66.61$\pm$2.56
Transformer Based Methods | Transformer [10] | 67.80$\pm$ 0.98 | 69.14$\pm$ 1.87 | 62.99$\pm$ 6.55 | 72.50$\pm$ 4.07
BrainNetTransformer [11] | 71.22$\pm$ 0.97 | 73.60$\pm$ 2.38 | 67.70$\pm$ 8.27 | 74.29$\pm$ 5.54
HyperGraphs (Ours) | HyperGraphGCN | 73.41$\pm$1.18 | 76.66$\pm$1.02 | 74.67$\pm$5.76 | 71.89$\pm$5.91
HyperGALE (Ours) | 75.34$\pm$ 0.47 | 77.03$\pm$ 1.85 | 76.39$\pm$ 4.84 | 73.91$\pm$ 5.58

Hypergraph Convolutions: Now, the incidence matrix $\mathcal{H}$ is passed to
a successive series of convolution layers. In order for the convolution
operation to proceed the propagation is done among the shared hyperedges and
higher confidence is assigned to the larger weights. Based on the hypergraph
structure and the associated weights, we define the hypergraph convolution
operation as,

$x_{t}^{(d+1)}=\sigma\left(\sum_{i=1}^{N}\sum_{j=1}^{K}\mathcal{H}_{tj}\mathcal{H}_{ij}W_{jj}x_{t}^{(d)}\Theta\right)$
(1)

Here $x_{t}^{(d)}$ is the embedding of the $t^{\text{th}}$ vertex in the
$d^{\text{th}}$ layer. The $\Theta$ is the weight matrix between two
successive layers. Now we can eventually formulate the equation (1) into a
matrix form as

$X^{(d+1)}=\sigma\left(\mathcal{H}W\mathcal{H}^{T}X^{(d)}\Theta\right)$ (2)

From the above equation the $X^{(d)}\in\mathbb{R}^{N\times F^{(d)}}$ is the
input of the $d^{\text{th}}$ layer and the $\Theta\in\mathbb{R}^{F^{(d)}\times
F^{(d+1)}}$. The equation (2) refers to a basic formulation that might lead to
numerical instabilities with deepening layers causing gradients to vanish.
Thus, to tackle this, one can adopt normalization approaches, such as
symmetric normalization or row-normalization. In our line of research, we
adopt the row-normalization approach as it provides directional propagation as
in the below equation

$X^{(d+1)}=\sigma\left(\mathcal{D}^{-1}\mathcal{H}W\mathcal{B}^{-1}\mathcal{H}^{T}X^{(d)}\Theta\right)$
(3)

Learnable Hyperedges and Gated Attention In the equation (3) the weight matrix
is not learnable. In such cases, the hypergraph will not be updated with
respect to the samples. So our method even tries to learn the hyperedge
weights matrix $W$ and we represent it as $\tilde{W}$ (refer to fig. 2).
Successively, this information is passed through a gated attention network as
formulated below,

$\displaystyle A=\sigma(\textsc{Mlp}(X^{(d+1)}));$ $\displaystyle
B=\text{Tanh}(\textsc{Mlp}(X^{(d+1)}));$ (4)
$\displaystyle\alpha=\textsc{MLP}(A\odot B);$ $\displaystyle Z=\alpha\odot
X^{(d+1)}$

Figure 3:  Illustration of the Gated Attention Module. GA module which learns
$\alpha$ iteratively (or numerically) from eq. (III-B), this $\alpha$ vector
is multiplied to get the final node features which is fed to the readout
layer. [33].

Where, $\odot$ determines the element-wise product and $\text{Tanh}(\cdot)$,
$\sigma(\cdot)$ are the nonlinear activation’s respectively. In the equation
(III-B), $Z$ obliges our model to allocate necessary attention that are
crucial to interpreting the exact incidence of ASD [33]. This guides the
hyperedge weight matrix to direct the importance to the essential nodes. Thus,
it eventually tracks essential regions of the brain that are actually
connected and injects the importance accordingly (refer to Fig. 3).

Readout Mechanisms The Readout methods produce an aggregation of node-level
features to form graph-level representations. The standard readouts such as,
$\textsc{max}(\cdot)$ and $\textsc{mean}(\cdot)$ underfit while learning on
brain hypergraphs. So we experimented with adaptive readouts like
SetTransformer and Janossy (permutation invariant) and MLP (permutation non-
invariant). Thus, we deploy _Neural Readouts_ where the hypothesis space of
the aggregated representations is not permutation invariant. Also similar to
Buterez _et al._ [34] we applied three variants of adaptive neural readouts,
SetTransformer, Janossy and Mlp readouts. We experimented with both
permutation variant and invariant that are adaptive and embedded with
permutation invariant learning with and without attention mechanisms [34].
After a detailed experiments we have seen Mlp readouts outperform many
adaptive and non-adaptive readouts. The ablations are detailed briefly in the
results section. Our chosen readout is,

$X^{\text{(Readout)}}=\textsc{Relu}(\textsc{Mlp}(Z))$ (5)

Finally, the total pipeline of our HyperGALE is described in the below
Algorithm.1

Input: Data

Output: $\hat{y}$

$FC\leftarrow\textsc{SchaeferParcellation}(\textsc{Data})$

for  _i $\leftarrow$ 1 to epochs_ do

$\;\;\text{Hypergraph Creation: }\mathcal{G}_{hyper}(V,\mathcal{H})$
$X^{(d+1)}\leftarrow\sigma\left(\mathcal{D}^{-1}\mathcal{H}\tilde{W}\mathcal{B}^{-1}\mathcal{H}^{T}X^{(d)}\Theta\right)$
$\;\;Z\leftarrow\alpha\odot X^{(d+1)}$

$\;\;X^{\text{(Readout)}}\leftarrow\textsc{Relu}(\textsc{Mlp}(Z))$

$\;\;\hat{y}\leftarrow\sigma(X^{\text{(Readout)}})$

return _$\hat{y}$_

Algorithm 1 HyperGALE

##  IV Results and Discussion

In this section, we compare our proposed method with methods that are employed
to provide significant performance for ASD classification. Later we determine
the qualitative performance of HyperGALE in determining ASD.

###  IV-A Methods for Comparison

All the baseline experiments have been conducted on the dataset we have
identified from ABIDE II and where necessary we have re-implemented the
algorithms for reporting the comparative results in Table I. First, to observe
the significance of statistical approaches the traditional machine learning
methods such as Support Vector Machines, Random Forests, and Gradient Boosting
Classifier (GBC) have been implemented. These models have been poorly
performed due to the high variance nature of the data (refer Table I).
Nevertheless, comparing among traditional methods, SVM gave better prediction
performance.

Next, moving to the architectures using deep learning without transformer, the
works [8, 9] have decent and better performance compared to traditional
methods. Although, these methods have succeeded to provide jump in performance
compared to that of traditonal methods these models have high standard
deviation and needs to be better in capturing the patterns that corresponds to
ASD.

Next, we move to latest transformer architectures [10, 11] have achieved
promising performance with less standard deviation in AUC and accuracy
metrics. Although, having self-attention layers and OCReadout mechansims,
these methods do not provide significant performance.

Now, we benchmark the graph-based models [13, 14, 15]. As explained earlier,
graph-based have actually performed well for ASD dataset [16][17] [18]. Graph-
based approaches on FC matrices are applied to understand interactions among
regions of interest, but current methods struggle with higher-order proximity
due to reliance on static dyadic edges.

Recent studies using Hypergraphs for ASD classification include Hypergraph
U-Net [19] and multi-view HGNN[20] which use multiple modalities and thus we
can’t employ these methods for our study. Finally, we employed the hypergraph
convolution methods named HyperGraphGCN [30] to compare with the hypergraph
baseline. Although these methods have the added advantage of acquiring higher-
order proximity, HyperGALE performs extensively on top of every model. In the
next subsection, we explain all the experiments in detail and justify
HyperGALE with extensive ablation studies.

TABLE II: The Table illustrates the performance hypergraph models in the
combinations of with and without GA and learned HE.

Methods | Performance Metrics
---|---
Learned HEs | GA | Accuracy | AUC | Sensitivity | Specificity
W/o | W/o | 73.41$\pm$1.18 | 76.66$\pm$1.02 | 74.67$\pm$5.76 | 71.89$\pm$5.91
W/o | With | 74.39$\pm$1.21 | 77.90$\pm$0.90 | 67.78$\pm$1.11 | 82.43$\pm$1.35
With | W/o | 73.74$\pm$1.39 | 77.15$\pm$3.80 | 68.10$\pm$7.20 | 81.15$\pm$6.56
With | With | 75.34$\pm$0.47 | 77.03$\pm$1.85 | 76.56$\pm$5.41 | 73.86$\pm$5.41

###  IV-B Performance Discussion

Table I showcases that HyperGALE achieves notable performance, effectively
justifying its design principles. The robustness of HyperGALE is confirmed in
both accuracy and AUC metrics. A two-sample t-test comparing HyperGALE with
HyperGraphGCN (a hypergraph-based baseline) reveals a significant difference
in accuracy (t-Statistic: 3.40, p-Value: 0.0094), highlighting HyperGALE’s
superior performance in this regard. When compared with BrainNetTransformer,
HyperGALE exhibits a highly significant improvement in accuracy (t-Statistic:
8.55, p-Value: 0.000027). Additionally, a significant enhancement is also
observed in AUC (t-Statistic: 2.54, p-Value: 0.0345). These results highlight
the efficacy of HyperGALE in handling functional connectivity (FC) patterns,
particularly in Autism Spectrum Disorder (ASD) pattern identification.

While GCN demonstrates promising sensitivity, it suffers from fluctuating
performance characterized by a high standard deviation. In terms of
specificity, BrainNetTransformer excels, yet HyperGALE achieves comparable
results. The superior performance of HyperGALE, especially in processing FC
patterns, emphasizes its capability in delineating nuanced ASD patterns.
Unlike transformer-based models that struggle with the utilization of FC as
tokens and the physical understanding of attention mechanisms, HyperGALE
adeptly processes these FC patterns, thereby capturing functional associations
crucial in ASD diagnostics. This approach not only bolsters the model’s
accuracy and AUC but also enriches the understanding of ASD through cutting-
edge machine learning techniques.

###  IV-C Generalizability

Our study employed a leave-one-site-out analysis across 16 different sites
(refer Fig. 4) in the ABIDE-II dataset to evaluate the HyperGALE model’s
generalizability. This approach rigorously tested the model’s performance in
diverse settings, revealing insights into its applicability across various
clinical environments. The analysis in Fig. 4 demonstrates varied performance
metrics - across sites.

The results underscore the challenge of generalizing models in multi-site
studies, highlighting the influence of site-specific factors such as data
acquisition protocols and demographic variations. Performance discrepancies
across sites suggest the need for further model enhancements, potentially
through site-adaptive methodologies.

Figure 4: Prediction performance across sites using a leave-one-site-out
strategy. The count of ASD and TD subjects from various sites are shown on
(left) and various performance metrics are shown for each site on (right).
Chance level performance is at 50%. Despite the challenges with site-specific
variations in the number of samples, our model is still able to demonstrate
creditable between-site generalization performance, comparable to the results
in Table I.

###  IV-D Gated Attention and Hyperedge Discussion

Table II demonstrates the impact of incorporating Learned Hyperedges (HE) and
the Gated Attention (GA) module in hypergraph models, assessed under various
configurations with and without these modules. One can observe that, with GA,
the model is performing well in AUC and specificity metrics even in the
absence of learnable HE. But, learnable HE and GA together actually obliges
the model to have higher accuracy score and sensitivity. The sensitivity is
helpful to predict the presence the ASD. Hence, both of these act as tandem
and pivotal in identifying key ROI embeddings, thereby focusing on the most
significant regions for predicting ASD. Furthermore, as GA highlights
significant ROIs, and learned HE assign varying weights to different networks
in brain parcellation, learning in this way helps the model to focus on
distinct ROIs to provide significant performance.

The significance of learned HE and the GA module extends beyond performance;
they are also crucial for enhancing the model’s interpretability in
identifying biomarkers, a topic we will explore in detail in Section V.

###  IV-E Ablation Study

#### Number of ROIs in each Hyperedge

In our hypergraph convolution approach, hyperedges are formed by grouping ROIs
based on connectivity, with the $k$-NN algorithm determining the number of
ROIs in each hyperedge. The choice of the hyperparameter $k$, which dictates
the ROI count per hyperedge, is crucial as it influences the overlap of
hyperedges, impacting the identification of active brain regions in ASD and TD
subjects [35]. A lower $k$ value might exclude vital information by limiting
the ROIs per hyperedge, whereas a higher $k$ risks adding non-significant
connections and noise. Our empirical analysis as seen in Fig. 5, determined 40
to be the optimal ROI count in a hyperedge, balancing comprehensive coverage
and minimal noise, crucial for accurate brain region representation in ASD and
TD subjects.

#### Number of Layers

Our analysis shows that a single-layer hypergraph convolution is optimal for
brain connectome analysis, as evidenced in Fig. 6. An increase in the number
of layers led to a decline in performance, a phenomenon consistent with over-
smoothing issues in graph convolutional networks (GCNs). This issue, as
discussed in [36, 37], involves feature homogenization across layers, reducing
the model’s ability to distinguish unique patterns.

TABLE III: The Table illustrates the performance of various readout mechanisms on our proposed HyperGALE model. Methods | Performance Metrics
---|---
Accuracy | AUC | Sensitivity | Specificity
MLP | 75.34 | 77.03 | 68.10 | 81.15
Set Transformer | 67.36 | 65.95 | 41.03 | 77.14
Janossy | 62.96 | 58.08 | 63.16 | 62.79
Max | 51.12 | 54.14 | 23.00 | 97.30
Mean | 48.17 | 50.36 | 10.00 | 94.59

The brain connectome’s structure, mainly consisting of ROIs one hop away from
each other, supports our use of a single-layer approach. This structure allows
for effective aggregation of essential information in the initial layer,
making additional layers redundant. While various strategies have been
proposed to address over-smoothing in deeper GNNs [38, 39], the architectural
specifics of our dataset make these approaches less relevant, as our single-
layer model adequately captures the necessary connectivity information.

#### Readout Mechanisms

Table III demonstrates that standard readouts like $\textsc{max}(\cdot)$ and
$\textsc{mean}(\cdot)$ are less effective for brain hypergraphs. Supporting
this, Buterez _et al._ [34] note that not all readouts suit every graph
structure, especially in fully connected graphs like brain hypergraphs where
permutation non-invariant readouts can be more effective.

The SetTransformer readout, used without positional encoding to preserve
permutation invariance, contrasts with the permutation-invariant Janossy
readout suggested by Murphy _et al._[40]. However, in our study, MLP readouts
provided superior embeddings for brain hypergraphs, aligning with findings by
Alon _et al._ [41] that demonstrate their effectiveness in similar contexts.
Hence, our framework utilized MLP readouts for optimal representation.

Figure 5: Changes in prediction accuracy with respect to change in number of
ROIs in a hyperedge denoted by k. The optimal number of ROIs was found at $k$
= 40. Figure 6: Changes in prediction accuracy with respect to change in the
number of hypergraph convolution layers. A single-layer hypergraph convolution
is found to be optimal.

##  V Interpretability Discussion

Our analysis of brain connectome interactions unveils distinct neural patterns
between ASD and TD subjects. Utilizing the Gated Attention module and learned
hyperedge weights, we identified pivotal ROIs and observed how their
interactions differ between ASD and TD individuals.

While Figs 7 (A) and (B) reveal a significant overlap in the top ROIs across
both ASD and TD groups, it is the distinct connections among regions of these
two groups that provide deeper insights and distinguish the two groups. In
other words, despite sharing common ROIs, the way these regions interact
within the neural network varies considerably between ASD and TD, reflecting
the unique neural connectivity characteristics of ASD.

We identified the top 50 ROIs and categorized them into broader functional
networks, each reflecting distinct aspects of the ASD classification. In the
following we highlight the differences observed in the functional networks
between ASD and TD. These include, Visual Processing Networks, which highlight
atypical visual processing patterns in ASD [42]. The Limbic System indicates
altered emotional and social processing through higher connectivity with the
Default Mode Network (DMN) [43]. Executive Function and Attention Networks
reveal disparities related to executive functioning and attention regulation
challenges in ASD; and the Default Mode Network, where atypical connectivity
patterns impact areas like theory of mind and self-awareness [43].

Figure 7:  The interpretation of important ROIs and their connectivity
obtained from GA module and learned hyperedge weights, respectively. (A) and
(B) highlight the top 20 most important ROIs that emerged in ASD and TD
subjects. (C) and (D) depict the _discriminative hyperedge connectivity
pattern_ of a _common_ ROI in ASD and TD using dominant hyperedge based on
$\tilde{W}$.

While Figs. 7 (A) and (B) indicate a substantial overlap in the top ROIs
across both ASD and TD groups, it is their distinctive connectivity patterns
that yield profound insights. This is particularly evident in the Limbic
System, comprising of regions like Limbic TempPole and Limbic OFC. In ASD, we
observe higher connectivity of these limbic regions with the DMN, as opposed
to the balanced connections seen in TD. This higher linkage in ASD, as
highlighted in Fig 7 (C) versus (D), is indicative of altered connectivity
patterns in areas crucial for emotional processing and social interaction
[43]. Similarly, in visual regions (RH Vis, LH Vis), we observe increased
connectivity to DMN regions in ASD subjects [42], while TD had connections
primarily within visual and attention networks.

These observations improve our understanding of ASD’s neural dynamics. The
Limbic System’s enhanced DMN connectivity in ASD might underlie social and
emotional processing challenges, while the varied connectivity in visual
networks points to broader sensory and cognitive implications. Overall, our
study extends beyond identifying shared ROIs, and discovering complex
variations in neural interactions.

##  VI Conclusion

In this paper, we introduced HyperGALE, a novel hypergraph convolutional
network enhanced with gated attention, specifically designed for ASD
classification. Our comprehensive evaluation on the ABIDE-II dataset
highlights HyperGALE’s outstanding performance, surpassing ML and graph-based
approaches. HyperGALE’s capability to effectively capture complex higher-order
graph intricacies and utilize learnable hyperedge weights has led to a better
understanding of brain network dynamics. HyperGALE identifies critical ROIs
for ASD and delineates how their connectivity patterns differ from those in TD
subjects. This, not only promises improvements in clinical detection but also
opens new avenues for research in neuroimaging and the development of targeted
intervention strategies.

## References

  * [1] M. J. Maenner et al., “Prevalence and characteristics of autism spectrum disorder among children aged 8 years—autism and developmental disabilities monitoring network, 11 sites, united states, 2020,” MMWR Surveillance Summaries, vol. 72, no. 2, p. 1, 2023.
  * [2] L. Q. Uddin, K. Supekar, and V. Menon, “Reconceptualizing functional brain connectivity in autism from a developmental perspective,” Frontiers in human neuroscience, vol. 7, p. 458, 2013.
  * [3] J. S. Nomi and L. Q. Uddin, “Developmental changes in large-scale network connectivity in autism,” NeuroImage: Clinical, vol. 7, pp. 732–741, 2015.
  * [4] S. M. Smith et al., “Correspondence of the brain’s functional architecture during activation and rest,” Proceedings of the national academy of sciences, vol. 106, no. 31, pp. 13040–13045, 2009.
  * [5] J. Bathelt, H. M. Geurts, and D. Borsboom, “More than the sum of its parts: Merging network psychometrics and network neuroscience with application in autism,” Network Neuroscience, vol. 6, no. 2, pp. 445–466, 2022.
  * [6] M. R. Arbabshirani, S. Plis, J. Sui, and V. D. Calhoun, “Single subject prediction of brain disorders in neuroimaging: Promises and pitfalls,” NeuroImage, vol. 145, pp. 137–165, 2017.
  * [7] Y. Chen, J. Yan, et al., “Adversarial learning based node-edge graph attention networks for autism spectrum disorder identification,” IEEE Transactions on Neural Networks and Learning Systems, p. 1–12, 2024.
  * [8] T. Eslami, V. Mirjalili, A. Fong, A. R. Laird, and F. Saeed, “Asd-diagnet: a hybrid learning approach for detection of autism spectrum disorder using fmri data,” Frontiers in neuroinformatics, vol. 13, p. 70, 2019.
  * [9] J. Kawahara et al., “Brainnetcnn: Convolutional neural networks for brain networks; towards predicting neurodevelopment,” NeuroImage, vol. 146, pp. 1038–1049, 2017.
  * [10] A. Vaswani et al., “Attention is all you need,” Advances in neural information processing systems, vol. 30, 2017.
  * [11] X. Kan, W. Dai, H. Cui, Z. Zhang, Y. Guo, and C. Yang, “Brain network transformer,” Advances in Neural Information Processing Systems, vol. 35, pp. 25586–25599, 2022.
  * [12] R.-A. Müller, N. Kleinhans, N. Kemmotsu, K. Pierce, and E. Courchesne, “Abnormal variability and distribution of functional maps in autism: an fmri study of visuomotor learning,” American Journal of Psychiatry, vol. 160, no. 10, pp. 1847–1862, 2003.
  * [13] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” arXiv preprint arXiv:1609.02907, 2016.
  * [14] P. Veličković et al., “Graph attention networks,” arXiv preprint arXiv:1710.10903, 2017.
  * [15] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” Advances in neural information processing systems, vol. 30, 2017.
  * [16] M. Cao, M. Yang, C. Qin, X. Zhu, Y. Chen, J. Wang, and T. Liu, “Using deepgcn to identify the autism spectrum disorder from multi-site resting-state data,” Biomedical Signal Processing and Control, vol. 70, p. 103015, 2021.
  * [17] A. Kazi, S. Shekarforoush, S. Arvind Krishna, H. Burwinkel, G. Vivar, K. Kortüm, S.-A. Ahmadi, S. Albarqouni, and N. Navab, “Inceptiongcn: receptive field aware graph convolutional network for disease prediction,” in Information Processing in Medical Imaging: 26th International Conference, IPMI 2019, Hong Kong, China, June 2–7, 2019, Proceedings 26, pp. 73–85, Springer, 2019.
  * [18] D. Yao, M. Liu, M. Wang, C. Lian, J. Wei, L. Sun, J. Sui, and D. Shen, “Triplet graph convolutional network for multi-scale analysis of functional connectivity using functional mri,” in Graph Learning in Medical Imaging: First International Workshop, GLMI 2019, Held in Conjunction with MICCAI 2019, Shenzhen, China, October 17, 2019, Proceedings 1, pp. 70–78, Springer, 2019.
  * [19] M. Lostar and I. Rekik, “Deep hypergraph u-net for brain graph embedding and classification,” arXiv preprint arXiv:2008.13118, 2020.
  * [20] M. Madine, I. Rekik, and N. Werghi, “Diagnosing autism using t1-w mri with multi-kernel learning and hypergraph neural network,” in 2020 IEEE International Conference on Image Processing (ICIP), pp. 438–442, IEEE, 2020.
  * [21] W. Shao, Y. Peng, C. Zu, M. Wang, D. Zhang, A. D. N. Initiative, et al., “Hypergraph based multi-task feature selection for multimodal classification of alzheimer’s disease,” Computerized Medical Imaging and Graphics, vol. 80, p. 101663, 2020.
  * [22] M. Liu, J. Zhang, P.-T. Yap, and D. Shen, “View-aligned hypergraph learning for alzheimer’s disease diagnosis with incomplete multi-modality data,” Medical image analysis, vol. 36, pp. 123–134, 2017.
  * [23] J. Ji, Y. Ren, and M. Lei, “Fc–hat: Hypergraph attention network for functional brain network classification,” Information Sciences, vol. 608, pp. 1301–1316, 2022.
  * [24] C. Craddock, S. Sikka, B. Cheung, R. Khanuja, S. S. Ghosh, C. Yan, Q. Li, D. Lurie, J. Vogelstein, R. Burns, S. Colcombe, M. Mennes, C. Kelly, A. Di Martino, F. X. Castellanos, and M. Milham, “Towards automated analysis of connectomes: The Configurable Pipeline for the Analysis of Connectomes (C-PAC),” Frontiers in Neuroinformatics, no. 42, 2013.
  * [25] B. B. Avants, N. J. Tustison, G. Song, P. A. Cook, A. Klein, and J. C. Gee, “A reproducible evaluation of ants similarity metric performance in brain image registration,” NeuroImage, vol. 54, no. 3, pp. 2033–2044, 2011.
  * [26] Y. Behzadi, K. Restom, J. Liau, and T. T. Liu, “A component based noise correction method (compcor) for bold and perfusion based fmri,” NeuroImage, vol. 37, no. 1, pp. 90–101, 2007.
  * [27] A. Schaefer, R. Kong, E. M. Gordon, T. O. Laumann, X.-N. Zuo, A. J. Holmes, S. B. Eickhoff, and B. T. T. Yeo, “Local-Global Parcellation of the Human Cerebral Cortex from Intrinsic Functional Connectivity MRI,” Cerebral Cortex, vol. 28, no. 9, pp. 3095–3114, 2017.
  * [28] O. Ledoit and M. Wolf, “A well-conditioned estimator for large-dimensional covariance matrices,” J. Multivar. Anal., vol. 88, p. 365, 2004.
  * [29] A. Abraham, F. Pedregosa, M. Eickenberg, P. Gervais, A. Mueller, J. Kossaifi, A. Gramfort, B. Thirion, and G. Varoquaux, “Machine learning for neuroimaging with scikit-learn,” Frontiers in neuroinformatics, vol. 8, 2014.
  * [30] Y. Feng, H. You, Z. Zhang, R. Ji, and Y. Gao, “Hypergraph neural networks,” in Proceedings of the AAAI conference on artificial intelligence, vol. 33, pp. 3558–3565, 2019.
  * [31] C. Cortes and V. Vapnik, “Support-vector networks,” Machine learning, vol. 20, pp. 273–297, 1995.
  * [32] L. Breiman, “Random forests,” Machine learning, vol. 45, pp. 5–32, 2001.
  * [33] M. Ilse, J. Tomczak, and M. Welling, “Attention-based deep multiple instance learning,” in International conference on machine learning, pp. 2127–2136, PMLR, 2018.
  * [34] D. Buterez, J. P. Janet, S. J. Kiddle, D. Oglic, and P. Liò, “Graph neural networks with adaptive readouts,” Advances in Neural Information Processing Systems, vol. 35, pp. 19746–19758, 2022.
  * [35] G. Lee, M. Choe, and K. Shin, “How do hyperedges overlap in real-world hypergraphs?-patterns, measures, and generators,” in Proceedings of the web conference 2021, pp. 3396–3407, 2021.
  * [36] Q. Li, Z. Han, and X.-M. Wu, “Deeper insights into graph convolutional networks for semi-supervised learning,” in Proceedings of the AAAI conference on artificial intelligence, vol. 32, 2018.
  * [37] D. Chen, Y. Lin, W. Li, P. Li, J. Zhou, and X. Sun, “Measuring and relieving the over-smoothing problem for graph neural networks from the topological view,” in Proceedings of the AAAI conference on artificial intelligence, vol. 34, pp. 3438–3445, 2020.
  * [38] M. Liu, H. Gao, and S. Ji, “Towards deeper graph neural networks,” in Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pp. 338–348, 2020.
  * [39] G. Li, M. Müller, G. Qian, I. C. D. Perez, A. Abualshour, A. K. Thabet, and B. Ghanem, “Deepgcns: Making gcns go as deep as cnns,” IEEE transactions on pattern analysis and machine intelligence, 2021.
  * [40] R. L. Murphy, B. Srinivasan, V. Rao, and B. Ribeiro, “Janossy pooling: Learning deep permutation-invariant functions for variable-size inputs,” arXiv preprint arXiv:1811.01900, 2018.
  * [41] U. Alon and E. Yahav, “On the bottleneck of graph neural networks and its practical implications,” International Conference on Learning Representations, 2021.
  * [42] M. V. Lombardo, L. Eyler, A. Moore, M. Datko, C. Carter Barnes, D. Cha, E. Courchesne, and K. Pierce, “Default mode-visual network hypoconnectivity in an autism subtype with pronounced social visual engagement difficulties,” Elife, vol. 8, p. e47427, 2019.
  * [43] E. D. Bigler, D. F. Tate, E. S. Neeley, L. J. Wolfson, M. J. Miller, S. A. Rice, H. Cleavinger, C. Anderson, H. Coon, S. Ozonoff, et al., “Temporal lobe, autism, and macrocephaly,” American Journal of Neuroradiology, vol. 24, no. 10, pp. 2066–2076, 2003.