# Supplementary Material for “Core-valence attosecond transient absorption
spectroscopy of polyatomic molecules”

Nikolay V. Golubev<EMAIL_ADDRESS>Laboratory of Theoretical Physical
Chemistry, Institut des Sciences et Ingénierie Chimiques, Ecole Polytechnique
Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Jiří Vaníček
Laboratory of Theoretical Physical Chemistry, Institut des Sciences et
Ingénierie Chimiques, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland  Alexander I. Kuleff  Theoretische Chemie, Universität
Heidelberg, Im Neuenheimer Feld 229, D-69120 Heidelberg, Germany  ELI-ALPS,
Wolfgang Sandner utca 3, H-6728 Szeged, Hungary

###### Abstract

In the Letter, we present an application of the attosecond transient
absorption spectroscopy (ATAS) to probe the electronic oscillations of a
positive charge created after outer-valence ionization of the propiolic acid
molecule, while accounting for all electronic and nuclear degrees of freedom.
Here, we present additional details: in Sec. I, we derive the expression for
the absorption cross-section in ATAS of polyatomic molecules. Section II
contains details of the procedure used for evaluating transition properties
between electronic states computed by high-level ab initio electronic
structure methods. In Secs. III and IV, we discuss the averaging of the
absorption cross-sections over rotational degrees of freedom and with respect
to relative contributions of ionic states to the initial wave packet,
respectively.

###### pacs:

Valid PACS appear here

##  I Transient absorption spectroscopy

The absorption cross-section can be computed as (Tannor, 2007, p. 411)

$\sigma(\omega)=-\frac{4\pi\omega}{c}\frac{\text{Im}[\tilde{P}^{*}(\omega)\tilde{E}(\omega)]}{|\tilde{E}(\omega)|^{2}},$
(1)

where $\omega$ is the photon energy, $c$ is the speed of light in vacuum,
$\tilde{E}(\omega)$ is the Fourier transform of the electric field
$\tilde{E}(\omega)=\mathcal{F}[E(t)]$, and $\tilde{P}(\omega)$ is the Fourier
transform of the polarization function

$P(t)=\langle\Psi(t)|\hat{\mu}|\Psi(t)\rangle,$ (2)

where $\hat{\mu}$ is the electric dipole operator, and $|\Psi(t)\rangle$ is
the time-dependent wavefunction. In order to calculate the field-induced
polarization response, we need to solve the time-dependent Schrödinger
equation

$i\frac{d|\Psi(t)\rangle}{dt}=\hat{H}(t)|\Psi(t)\rangle,$ (3)

with the Hamilton operator

$\hat{H}(t)=\hat{H}_{0}-\hat{\mu}E(t),$ (4)

describing a system under the influence of an external field. Here, the scalar
dipole operator $\hat{\mu}=\hat{\vec{\mu}}\cdot\vec{\epsilon}$ denotes the
projection of the dipole vector operator $\hat{\vec{\mu}}$ onto the unit
polarization vector $\vec{\epsilon}$ of an external field
$\vec{E}(t)=\vec{\epsilon}\ E(t)$.

To evaluate the polarization $P(t)$, we employ the time-dependent perturbation
theory, in which the state at time $t$ is computed as

$|\Psi(t)\rangle=\hat{U}(t)|\Psi(0)\rangle+i\int_{-\infty}^{t}dt^{\prime}\hat{U}(t-t^{\prime})\hat{\mu}E(t^{\prime})\hat{U}(t^{\prime})|\Psi(0)\rangle+\cdots,$
(5)

where $\hat{U}(t)$ is the field-free evolution operator for the Hamiltonian
$\hat{H}_{0}$ and $|\Psi(0)\rangle$ is the initial state of the system at time
zero.

Using Eq. (5), we can expand the polarization $P(t)$ in powers of the electric
field as

$P(t)=P^{(0)}(t)+P^{(1)}(t)+\cdots,$ (6)

where the zeroth-order term $P^{(0)}(t)$ describes the harmonic emission,
resulting from oscillations of the dipole moment in time (see, e.g., Ref.
Kuleff and Cederbaum (2011)), and

$P^{(1)}(t)=i\langle\Psi(0)|\hat{U}^{\dagger}(t)\hat{\mu}\int_{-\infty}^{t}dt^{\prime}\hat{U}(t-t^{\prime})\hat{\mu}E(t^{\prime})\hat{U}(t^{\prime})|\Psi(0)\rangle+\text{c.c.},$
(7)

is the linear polarization response with respect to the applied electric
field. Importantly, in an isotropic medium even-order polarization terms (i.e.
$P^{(0)}(t)$, $P^{(2)}(t)$, etc.) do not survive orientational averaging which
will be employed later and thus can be ignored Tannor (2007). The higher-order
terms of the polarization expansion, Eq. (6), are responsible for various
processes that enter in a broad variety of non-linear spectroscopy techniques
Mukamel (1995), which are beyond the scope of the present study and thus will
be neglected. Henceforth, hereafter we will concentrate on the linear
polarization response function, Eq. (7), which describes absorption and
stimulated emission driven by the external electric field $E(t)$.

Let us assume for the moment that the eigenvalues $\varepsilon_{k}$ and
eigenvectors $|\Psi_{k}\rangle$ of the operator $\hat{H}_{0}$ are known. Since

$\hat{H}_{0}|\Psi_{k}\rangle=\varepsilon_{k}|\Psi_{k}\rangle,$ (8)

the action of the evolution operator $\hat{U}(t)$ on the initial state
$|\Psi(0)\rangle$ can be represented as

$\hat{U}(t)|\Psi(0)\rangle=\sum_{k}c_{k}e^{-i\varepsilon_{k}t}|\Psi_{k}\rangle,$
(9)

where $c_{k}=\langle\Psi_{k}|\Psi(0)\rangle$ are the expansion coefficients of
the initial state in a basis of $\hat{H}_{0}$ eigenstates. Using the explicit
form of the evolution operator acting on the initial state in Eq. (9), we can
write the first-order polarization correction as

$P^{(1)}(t)=i\int_{-\infty}^{t}dt^{\prime}E(t^{\prime})\sum_{k}\sum_{j}c^{*}_{k}c_{j}e^{i\varepsilon_{k}t}e^{-i\varepsilon_{j}t^{\prime}}\langle\Psi_{k}|\hat{\mu}\hat{U}(t-t^{\prime})\hat{\mu}|\Psi_{j}\rangle+\text{c.c.}$
(10)

To evaluate the matrix elements
$\langle\Psi_{k}|\hat{\mu}\hat{U}(t-t^{\prime})\hat{\mu}|\Psi_{j}\rangle$, we
insert the resolution of identity
$\sum_{f}|\Psi_{f}\rangle\langle\Psi_{f}|\equiv 1$ between
$\hat{U}(t-t^{\prime})$ and $\hat{\mu}$ operators, thus obtaining

$P^{(1)}(t)=i\sum_{k}\sum_{j}c^{*}_{k}c_{j}\sum_{f}\langle\Psi_{k}|\hat{\mu}|\Psi_{f}\rangle\langle\Psi_{f}|\hat{\mu}|\Psi_{j}\rangle\underbrace{\int_{-\infty}^{t}dt^{\prime}E(t^{\prime})e^{i\varepsilon_{k}t}e^{-i\varepsilon_{j}t^{\prime}}e^{-i\varepsilon_{f}(t-t^{\prime})}}_{I(t)}+\text{c.c.}$
(11)

By reorganizing $\varepsilon$ energies on the right-hand side of Eq. (11), the
time integral $I(t)$ can be represented in the form of a convolution

$I(t)=\int_{-\infty}^{t}dt^{\prime}E(t^{\prime})e^{i(\varepsilon_{k}-\varepsilon_{j})t^{\prime}}e^{-i(\varepsilon_{f}-\varepsilon_{k})(t-t^{\prime})}=\int_{-\infty}^{\infty}dt^{\prime}E(t^{\prime})e^{i(\varepsilon_{k}-\varepsilon_{j})t^{\prime}}e^{-i(\varepsilon_{f}-\varepsilon_{k})(t-t^{\prime})}\theta(t-t^{\prime}),$
(12)

where $\theta(t-t^{\prime})$ is the Heaviside step function. The convolution
of two functions in the time domain gives a simple product in the frequency
domain (Fourier convolution theorem)

$\tilde{I}(\omega)=\sqrt{2\pi}\mathcal{F}[E(t)e^{i(\varepsilon_{k}-\varepsilon_{j})t}]\mathcal{F}[e^{-i(\varepsilon_{f}-\varepsilon_{k})t}\theta(t)],$
(13)

where both terms in the product can be evaluated analytically:

$\mathcal{F}[E(t)e^{i(\varepsilon_{k}-\varepsilon_{j})t}]=\tilde{E}(\omega-\epsilon_{k}+\epsilon_{j}),$
(14)

and

$\mathcal{F}[e^{-i(\varepsilon_{f}-\varepsilon_{k})t}\theta(t)]=\frac{1}{\sqrt{2\pi}}\frac{-i}{\epsilon_{f}-\epsilon_{k}-\omega}.$
(15)

Importantly, the energy difference $\epsilon_{f}-\epsilon_{k}$ between states
coupled by the electric field has to be taken complex with a negative
imaginary part in order to converge the Fourier transform in Eq. (15). From a
physical point of view, complex values of energy levels can be associated with
a finite lifetime of final states which will result in the broadening of the
absorption lines in the spectrum.

Approximating the laser field by a delta function centered at $t=\tau$, i.e.,
$E(t)=E_{0}\delta(t-\tau)$, the Fourier transform $\tilde{E}(\omega)$ of the
electric field yields

$\tilde{E}(\omega)=\frac{1}{\sqrt{2\pi}}E_{0}e^{i\omega\tau},$ (16)

while the Fourier transform in Eq. (14) becomes

$\tilde{E}(\omega-\epsilon_{k}+\epsilon_{j})=\frac{1}{\sqrt{2\pi}}E_{0}e^{i(\omega+\epsilon_{k}-\epsilon_{j})\tau}.$
(17)

Substituting the above presented relations into Eq. (1), we obtain the
following expression for the absorption cross-section (see also Ref. Santra
_et al._ (2011))

$\sigma(\omega,\tau)=\frac{4\pi\omega}{c}\text{Im}\sum_{k}\sum_{j}c^{*}_{k}e^{i\varepsilon_{k}\tau}c_{j}e^{-i\varepsilon_{j}\tau}\sum_{f}\langle\Psi_{k}|\hat{\mu}|\Psi_{f}\rangle\langle\Psi_{f}|\hat{\mu}|\Psi_{j}\rangle\left(\frac{1}{\varepsilon_{f}-\varepsilon_{k}-\omega}+\frac{1}{\varepsilon_{f}^{*}-\varepsilon_{j}+\omega}\right).$
(18)

So far, we did not specify the physical nature of the Hamiltonian
$\hat{H}_{0}$ and the corresponding eigenstates and eigenenergies, Eq. (8).
Let us consider a particular case, where the operator $\hat{H}_{0}$ is the
full molecular Hamiltonian

$\hat{H}_{0}=\hat{T}_{n}+\hat{H}_{e},$ (19)

where $\hat{T}_{n}$ is the kinetic energy operator of the nuclei and
$\hat{H}_{e}$ denotes the electronic Hamiltonian. Using the well-known
Born–Oppenheimer formalism, i.e. solving the eigenvalue equation
$\hat{H}_{e}\Phi_{I}(\mathbf{r},\mathbf{R})=E_{I}(\mathbf{R})\Phi_{I}(\mathbf{r},\mathbf{R})$,
we can replace the summation over molecular states in Eq. (18) with a
summation over the electronic states:

$\sum_{j}c_{j}e^{-i\varepsilon_{j}\tau}|\Psi_{j}\rangle=\sum_{J}\chi_{J}(\mathbf{R},\tau)|\Phi_{J}\rangle,$
(20)

where $\chi_{J}(\mathbf{R},\tau)$ are nuclear wave packets propagating on the
corresponding potential energy surfaces $E_{J}(\mathbf{R})$. Furthermore, we
assume that vibrationally-resolved absorption lines in Eq. (18) can be
approximated with good accuracy by the averaged transitions between initial
$\epsilon_{j}\approx E_{J}$ and final
$\epsilon_{f}\approx\tilde{E}_{F}=E_{F}-i\frac{\Gamma}{2}$ electronic energy
levels, energetically broadened by the parameter $\Gamma$. The particular
value of $\Gamma$ has to be taken such that it accounts for the core-hole
lifetimes ($\sim$160 meV for the O${}_{\text{1s}}$ Sankari  _et al._ (2003)
and $\sim$100 meV for C${}_{\text{1s}}$ Carroll  _et al._ (2000)) and the
experimental resolution in ATAS (typically on the order of 0.2-0.3 eV at the
X-ray energy range employed in the present work (see, e.g., Ref. Geneaux  _et
al._ (2019))). The final form of the cross-section for electronically resolved
linear absorption spectroscopy reads

$\begin{split}\sigma(\omega,\tau)=\frac{4\pi\omega}{c}&\text{Im}\sum_{I}\sum_{J}\langle\chi_{I}(\mathbf{R},\tau)|\chi_{J}(\mathbf{R},\tau)\rangle_{\mathbf{R}}\sum_{F}\langle\Phi_{I}|\hat{\mu}|\Phi_{F}\rangle\langle\Phi_{F}|\hat{\mu}|\Phi_{J}\rangle\\\
&\times\left(\frac{1}{\tilde{E}_{F}-E_{I}-\omega}+\frac{1}{\tilde{E}^{*}_{F}-E_{J}+\omega}\right),\end{split}$
(21)

where the dynamics of initially populated states is driven by the electronic
coherence terms
$\langle\chi_{I}(\mathbf{R},\tau)|\chi_{J}(\mathbf{R},\tau)\rangle_{\mathbf{R}}$,
and the absorption lines occur between electronic states coupled by electric
dipole transitions.

##  II Transition properties between excited states

To calculate transition properties between states computed at different levels
of ab initio electronic structure theory, we employed a general technique
exploiting the configuration-interaction (CI)-like structure of the electronic
wavefunctions. The starting point is a CI expansion of an $N$-electron
wavefunction in terms of various excitations taking place from the reference
state $|\Psi_{0}\rangle$

$|\Phi_{I}\rangle=\sum_{p}a^{I}_{p}\hat{C}_{p}|\Psi_{0}\rangle,$ (22)

where $a^{I}_{p}$ are, in general, complex expansion coefficients, and
$\hat{C}_{p}$ is a formal representation of an operator describing excitations
in terms of creation $\hat{c}^{\dagger}$ and annihilation $\hat{c}$ operators
acting on the reference state. The representation of the wavefunction in the
form of Eq. (22) is general and appears in many methods of electronic
structure theory such as CI of various types Helgaker  _et al._ (2000),
equation-of-motion coupled-clusters (EOM-CC) family of methods Krylov (2008),
and algebraic diagrammatic construction (ADC) scheme Schirmer (2018), to name
a few.

A particular excitation $p$ in Eq. (22) represents a manipulation with
electrons occupying certain orbitals of the system which can be conveniently
described in terms of Slater determinants. Mathematically, the Slater
determinant of the reference state $|\Psi_{0}\rangle$ is defined as an
antisymmetrized product of orthonormal one-electron molecular spin-orbitals
$\\{\varphi_{1},\varphi_{2},...,\varphi_{N}\\}$

$\displaystyle|\Psi_{0}\rangle=\frac{1}{\sqrt{N!}}\begin{vmatrix}\varphi_{1}(\mathbf{x}_{1})&\varphi_{2}(\mathbf{x}_{1})&\cdots&\varphi_{N}(\mathbf{x}_{1})\\\
\varphi_{1}(\mathbf{x}_{2})&\varphi_{2}(\mathbf{x}_{2})&\cdots&\varphi_{N}(\mathbf{x}_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\varphi_{1}(\mathbf{x}_{N})&\varphi_{2}(\mathbf{x}_{N})&\cdots&\varphi_{N}(\mathbf{x}_{N})\end{vmatrix}\equiv|\varphi_{1},\varphi_{2},...,\varphi_{N}\rangle,$
(23)

where $\mathbf{x}_{i}$ denotes the spin and spatial coordinates of the $i$-th
electron. Let us consider an excitation
$\hat{C}_{p}=\hat{c}^{\dagger}_{a}\hat{c}_{k}$ generating one-hole–one-
particle (1$h$-1$p$) configurations: an electron is removed from an occupied
spin-orbital $\varphi_{k}$ and an electron is added to a virtual spin-orbital
$\varphi_{a}$

$|\Psi_{p}\rangle=\hat{C}_{p}|\Psi_{0}\rangle=\hat{c}^{\dagger}_{a}\hat{c}_{k}|\varphi_{1}...\varphi_{k}...\varphi_{N}\rangle=|\varphi_{1}...\varphi_{a}...\varphi_{N}\rangle,$
(24)

where in the singly excited Slater determinant $|\Psi_{p}\rangle$ the spin-
orbital $\varphi_{k}$ is replaced by another spin-orbital $\varphi_{a}$.
Similarly, any kind of electronic excitation $p$ can be represented by adding
or removing spin-orbitals from the reference Slater determinant, Eq. (23).

Representing the wavefunctions from Eq. (22) in terms of Slater determinants
allows us to express the matrix elements of a one-electron operator $\hat{F}$
between excited states $I$ and $J$ as

$\langle\Phi_{I}|\hat{F}|\Phi_{J}\rangle=\sum_{p}\sum_{q}(a^{I}_{p})^{*}a^{J}_{q}\langle\Psi_{p}|\hat{F}|\Psi_{q}\rangle,$
(25)

being proportional to the sum of all matrix elements
$\langle\Psi_{p}|\hat{F}|\Psi_{q}\rangle$ between the Slater determinants $p$
and $q$. These matrix elements can be expressed in terms of integrals in the
spin-orbital space using Slater–Condon rules Slater (1929); Condon (1930) (see
also Ref. Helgaker  _et al._ (2000) for a general review). For a one-electron
operator, these rules are:

$\langle\varphi_{1},\varphi_{2},...,\varphi_{N}|\hat{F}|\varphi_{1},\varphi_{2},...,\varphi_{N}\rangle=\sum_{i=1}^{N}\langle\varphi_{i}|\hat{F}|\varphi_{i}\rangle,$
(26)

i.e., the expectation value of the operator $\hat{F}$ between identical Slater
determinants is obtained by summing up matrix elements
$\langle\varphi_{i}|\hat{F}|\varphi_{i}\rangle$ between spin-orbitals;

$\langle\varphi_{1}...\varphi_{k}...\varphi_{N}|\hat{F}|\varphi_{1}...\varphi_{a}...\varphi_{N}\rangle=\langle\varphi_{k}|\hat{F}|\varphi_{a}\rangle,$
(27)

i.e. the matrix element between two Slater determinants differing by a single
orbital is equal to the integral between corresponding spin-orbitals; and

$\langle\varphi_{1}...\varphi_{k}...\varphi_{j}...\varphi_{N}|\hat{F}|\varphi_{1}...\varphi_{a}...\varphi_{b}...\varphi_{N}\rangle=0,$
(28)

when the Slater determinants differ by two or more spin-orbitals. Therefore,
computation of matrix elements of a one-electron operator between excited
states boils down to evaluating transition integrals between spin-orbitals of
the system.

Importantly, a particular form of the spin-orbital basis depends on the nature
of the reference state $|\Psi_{0}\rangle$ from which the excitations are
taking place. A simplest case is a standard CI scheme where the reference
state is the Hartree–Fock wavefunction, i.e. the spin-orbitals are obtained by
the diagonalization of the Fock matrix, and the state is a single Slater
determinant, Eq. (23). In more sophisticated methods such as EOM-CC or ADC,
the reference state is correlated and thus additional transformations of the
spin-orbital basis are required. In the case of ADC scheme used in the present
work, the reference state is constructed by applying the Rayleigh–Schrödinger
perturbation theory to the Hartree–Fock wavefunction and supposing a standard
Møller–Plesset partitioning of the Hamiltonian Schirmer (2018). Since the
resulting expressions, Eqs. (25) to (28), for evaluating the one-particle
operator $\hat{F}$ do not contain the spin-orbitals themselves but only the
matrix elements of the corresponding operator, the application of the
perturbation theory to the reference state is equivalent to the expansion of
the operator $\hat{F}$ in a perturbation series

$\mathbf{F}=\mathbf{F}^{(0)}+\mathbf{F}^{(1)}+\mathbf{F}^{(2)}+\cdots,$ (29)

where $\mathbf{F}$ denotes representation of the operator $\hat{F}$ in a basis
of spin-orbitals. The derivation of the explicit expressions for computing the
corresponding contributions $\hat{F}^{(i)}$ is beyond the scope of the present
work. However, we would like to note that the required matrix elements
$F_{kl}$ of the one-particle operator up to the second order of perturbations
can be straightforwardly obtained by combining the expressions for the
hole/hole and particle/hole terms given in appendices of Refs. Schirmer and
Trofimov (2004); Trofimov and Schirmer (2005).

##  III Orientational averaging of the transient absorption spectrum

Expression (21) is derived for a particular orientation of a molecule with
respect to the applied field. In the situation when the molecules are oriented
randomly in space, the cross-section can be obtained by averaging over all
possible orientations with respect to the laser polarization (see also Ref.
Begušić  _et al._ (2018) and references therein):

$\begin{split}\sigma(\omega,\tau)=&\int_{0}^{2\pi}\int_{0}^{\pi}\vec{\epsilon}(\theta,\varphi)^{T}\cdot\bm{\sigma}_{\vec{\mu},\vec{\mu}}(\omega,\tau)\cdot\vec{\epsilon}(\theta,\varphi)d\theta
d\varphi\\\
=&\frac{1}{3}\left(\sigma_{\mu_{x},\mu_{x}}(\omega,\tau)+\sigma_{\mu_{y},\mu_{y}}(\omega,\tau)+\sigma_{\mu_{z},\mu_{z}}(\omega,\tau)\right),\end{split}$
(30)

where
$\vec{\epsilon}(\theta,\varphi)^{T}=\left(\sin{(\theta)}\cos{(\varphi)},\
\sin{(\theta)}\sin{(\varphi)},\ \cos{(\theta)}\right)$ is the radial unit
vector, $\bm{\sigma}_{\vec{\mu},\vec{\mu}}$ is a tensor obtained by combining
different orthogonal projections of the dipole operator $\hat{\vec{\mu}}$ in
the corresponding transitions $\langle\Phi_{I}|\hat{\mu}|\Phi_{F}\rangle$ and
$\langle\Phi_{F}|\hat{\mu}|\Phi_{J}\rangle$ present in Eq. (21), and
$\sigma_{\mu_{x},\mu_{x}}(\omega,\tau)$,
$\sigma_{\mu_{y},\mu_{y}}(\omega,\tau)$,
$\sigma_{\mu_{z},\mu_{z}}(\omega,\tau)$ are absorption cross-sections computed
for three mutually orthogonal orientations of the molecule.

(a)

(b) (a) Absorption in $x$ direction. (c) (b) Absorption in $y$ direction.

Figure 1: Time-resolved absorption cross-section as a function of the photon
energy and time delay along two orthogonal directions of propiolic acid
molecule: $x$ (left) and $y$ (right). Top panels: Electronic coherence
measured by the time-dependent overlap
$\langle\chi_{1}(\mathbf{R},\tau)|\chi_{3}(\mathbf{R},\tau)\rangle_{\mathbf{R}}$
of the nuclear wave packets propagated in the first and third cationic states
of propiolic acid after the removal of an electron from the HOMO. Middle
panel: The absorption cross-section plotted for the energy window
corresponding to transitions between initially populated valence ionic states
and the core states resulting from ionization out of $1s$ orbitals of carbon
atoms forming the triple bond. Bottom panel: The absorption cross-section
plotted for the energy range covering transitions to core ionic states of
oxygen atoms. The ionic states involved in the corresponding transitions are
also shown.

The absorption cross-sections for different orientations of the propiolic acid
molecule with respect to the polarization of the applied field are shown in
Fig. 1. Due to the planar geometry of the molecule and the fact that the
initial superposition of the electronic states is prepared by ionizing an
electron from the HOMO belonging to A′ symmetry (see Fig. 1 of the main text),
the electric dipole transitions oriented perpendicularly to the molecular
plane vanish leading thus to complete transmission of the electric field
polarized in this direction. In contrast, a laser pulse polarized in the
molecular plane is absorbed in such a way that the maximum of the absorption
from initially created superposition of valence states to core ionic states
coincides with the localization of the charge density in the vicinity of a
particular atom. As one can see, the absorption on the triple-bond carbon
atoms takes place almost exclusively along $x$ direction of the molecule while
only a weak signal appears along $y$ direction. Indeed, the orbitals depicted
in Fig. 1 of the main text suggest that the positive and negative
contributions of the HOMO (shown in red and blue colors in Fig. 1 of the main
text) in the vicinity of the triple-bond along $y$ direction of the molecule
cancel each other out leading thus to nearly zero dipole moments connecting
the HOMO with the core orbitals of the carbon atoms in this direction. It is
also seen from the same figure that the HOMO-2 is mostly localized in the
vicinity of the carboxyl group which results in the weak, yet distinguishable,
dipole couplings with the core orbitals of the carbon atoms along both
directions. The stronger dipole transitions to the carbon orbitals along $y$
direction from the HOMO-2 in comparison with those from the HOMO is reflected
in the delay of the maximum of the absorption signal depicted in the right
panel of Fig. 1 of SM. Following similar logic, we can see from Fig. 1 of the
main text that the dipole transitions to the oxygen core orbitals are expected
to be stronger from the HOMO-2 in comparison with those from the HOMO. Indeed,
it is seen from Fig. 1 of SM that the maximum of the absorption involving
oxygen core orbitals is reached when the charge initially created in the HOMO
migrates to the HOMO-2. The averaged absorption cross-section obtained through
Eq. (30) is shown in Fig. 2 of the main text.

##  IV Averaging of the transient absorption spectrum with respect to the
relative contributions of ionic states to the initial wave packet

The sudden ionization limit employed in the present work assumes that the
ionization is performed with a short high-energy laser pulse. Due to a large
bandwidth of such a pulse, one can expect that the initial wave packet will be
prepared by a coherent superposition of multiple ionic states covered by the
applied pulse. Importantly, the large energy gap between the four valence
ionic states of the propiolic acid shown in Fig. 1 of the main text and the
ionic states lying higher in energy makes it possible to disentangle the high-
energy contributions from the transient absorption spectrum presented in this
work. Indeed, the electronic coherences between the lowest four cationic
states and the higher states will have oscillation periods well below a
femtosecond and can be removed from the spectrum, for example, by Fourier
analysis of the ATAS signal. Furthermore, the second and the fourth ionic
states shown in Fig. 1 of the main text are composed almost exclusively from
the one-hole configurations that overlap neither with each other, nor with the
configurations building the first and third cationic states of interest. The
lack of spatial overlap between the states in the created wave packet leads to
trivial phase dynamics, while observable properties of the system, such as the
charge density, will remain stationary in time. Therefore, the absorption
lines corresponding to transitions from the second and fourth ionic states to
the core ionic states will not oscillate and thus can be filtered out from the
spectrum.

The instability of laser pulse parameters, as well as the random orientation
of molecules in space might lead to variations of the populations of the ionic
states in the initial superposition. The absorption cross-sections for
different relative contributions of ionic states to the initial wave packet
created after ionization of the propiolic acid are shown in Fig. 2. Left and
middle panels of Fig. 2 illustrate the two limiting situations when the
initial state is prepared exclusively from the first or the third cationic
state, respectively. It is seen that in both cases the absorption cross-
sections remain stationary in time, whereby only the lines corresponding to
transitions from the first or the third valence state to the core ones are
present. Right panel of Fig. 2 demonstrates the averaged absorption cross-
section obtained by varying the relative contributions of the first and the
third ionic states to the initial wave packet, i.e.,
$\chi_{1}(\mathbf{R},0)=a\chi_{0}(\mathbf{R})$ and
$\chi_{3}(\mathbf{R},0)=\sqrt{1-a^{2}}\chi_{0}(\mathbf{R})$, where
$\chi_{0}(\mathbf{R})$ is the neutral ground state of the molecule. The
resulting spectrum is obtained by arithmetical averaging of ATAS signals
computed for eleven values of parameter $a$ taken between 0 and 1. As one can
see, the ATAS signal survives the averaging and thus provides sufficient
information to trace temporal and spatial aspects of the underlying electron
dynamics.

(a) (a) Only first cationic state is populated. (b) (b) Only third cationic
state is populated. (c) (c) Averaging over various populations of the first
and the third cationic states.

Figure 2: Time-resolved absorption cross-section as a function of the photon
energy and time delay obtained for different relative contributions of ionic
states to the initial wave packet created after ionization of the propiolic
acid. Left panel: Only the first cationic state is populated. Middle panel:
Only the third cationic state is populated. Right panel: The averaged spectrum
obtained by varying the relative contributions of the first and the third
ionic states in the initial wave packet. Top and bottom panels in each case
represent energy windows corresponding to transitions between initially
populated valence ionic states and the core states resulting from ionization
out of $1s$ orbitals of carbon and oxygen atoms, respectively.

## References

  * Tannor (2007) D. J. Tannor,  _Introduction to Quantum Mechanics: A Time-Dependent Perspective_ (University Science Books, Sausalito, 2007).
  * Kuleff and Cederbaum (2011) A. I. Kuleff and L. S. Cederbaum, Physical Review Letters 106, 053001 (2011).
  * Mukamel (1995) S. Mukamel,  _Principles of Nonlinear Optical Spectroscopy_ , Oxford Series in Optical and Imaging Sciences No. 6 (Oxford University Press, New York, 1995).
  * Santra  _et al._ (2011) R. Santra, V. S. Yakovlev, T. Pfeifer, and Z.-H. Loh, Physical Review A 83, 033405 (2011).
  * Sankari  _et al._ (2003) R. Sankari, M. Ehara, H. Nakatsuji, Y. Senba, K. Hosokawa, H. Yoshida, A. De Fanis, Y. Tamenori, S. Aksela, and K. Ueda, Chemical Physics Letters 380, 647 (2003).
  * Carroll  _et al._ (2000) T. X. Carroll, J. Hahne, T. D. Thomas, L. J. Sæthre, N. Berrah, J. Bozek, and E. Kukk, Physical Review A 61, 042503 (2000).
  * Geneaux  _et al._ (2019) R. Geneaux, H. J. B. Marroux, A. Guggenmos, D. M. Neumark, and S. R. Leone, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377, 20170463 (2019).
  * Helgaker  _et al._ (2000) T. Helgaker, P. Jørgensen, and J. Olsen,  _Molecular Electronic-Structure Theory_ (Wiley, New York, 2000).
  * Krylov (2008) A. I. Krylov, Annual Review of Physical Chemistry 59, 433 (2008).
  * Schirmer (2018) J. Schirmer,  _Many-Body Methods for Atoms, Molecules and Clusters_, Lecture Notes in Chemistry, Vol. 94 (Springer International Publishing, Cham, 2018).
  * Slater (1929) J. C. Slater, Physical Review 34, 1293 (1929).
  * Condon (1930) E. U. Condon, Physical Review 36, 1121 (1930).
  * Schirmer and Trofimov (2004) J. Schirmer and A. B. Trofimov, The Journal of Chemical Physics 120, 11449 (2004).
  * Trofimov and Schirmer (2005) A. B. Trofimov and J. Schirmer, The Journal of Chemical Physics 123, 144115 (2005).
  * Begušić  _et al._ (2018) T. Begušić, A. Patoz, M. Šulc, and J. Vaníček, Chemical Physics 515, 152 (2018).


Normative modeling has recently been introduced as a promising approach for modeling variation of neuroimaging measures across individuals in order to derive biomarkers of psychiatric disorders. Current implementations rely on Gaussian process regression, which provides coherent estimates of uncertainty needed for the method but also suffers from drawbacks including poor scaling to large datasets and a reliance on fixed parametric kernels. In this paper, we propose a deep normative modeling framework based on neural processes (NPs) to solve these problems. To achieve this, we define a stochastic process formulation for mixed-effect models and show how NPs can be adopted for spatially structured mixed-effect modeling of neuroimaging data. This enables us to learn optimal feature representations and covariance structure for the random-effect and noise via global latent variables. In this scheme, predictive uncertainty can be approximated by sampling from the distribution of these global latent variables. On a publicly available clinical fMRI dataset, we compare the novelty detection performance of multivariate normative models estimated by the proposed NP approach to a baseline multi-task Gaussian process regression approach and show substantial improvements for certain diagnostic problems.

Neural Processes, Mixed-Effect Modeling, Deep Learning, Neuroimaging.

§ INTRODUCTION

Recently, there has been great interest in applying machine learning to neuroimaging in order to find structural or functional biomarkers for brain disorders [Bzdok and Meyer-Lindenberg, 2018]. Such biomarkers can potentially be used for diagnosis or predicting treatment outcome in the spirit of precision medicine [Mirnezami et al., 2012]. In psychiatry, this is very challenging because clinical groups are highly heterogeneous in terms of symptoms and underlying biology  [Kapur et al., 2012]. However, most common analysis approaches ignore such heterogeneity and, in a case-control setting consider groups as distinct entities [Foulkes and Blakemore, 2018], where subjects are simply labeled as `patients' or `controls'. Supervised machine learning methods have been widely used in such settings but their accuracy is limited by the heterogeneity within each disorder [Wolfers et al., 2015].

Normative modeling [Marquand et al., 2016] is an emerging approach to address this challenge that has shown significant promise in multiple clinical settings [Wolfers et al., 2018, Zabihi et al., 2018, Wolfers et al., 2019]. Normative modeling involves estimating variation across the population in terms of mappings between clinically relevant covariates (e.g., age, cognitive scores) and biology (e.g., neuroimages). This is analogous to the use of `growth charts' in pediatric medicine to map variation in height or weight as a function of age. Currently, this is implemented using probabilistic regression methods that provide estimates of predictive uncertainty which map variation across the population. Deviations from the resulting normative model can then be interpreted as subject-specific biomarkers for brain disorders. For example, these can be used in a novelty detection setting for predicting diagnosis in an unsupervised fashion [Kia and Marquand, 2018, Kia et al., 2018].

Accurate quantification of uncertainty is crucial for normative modeling. In the original framework [Marquand et al., 2016], Gaussian process regression [Williams and Rasmussen, 1996] (GPR) was the central tool used to regress neuroimaging measures from clinical covariates. GPR is appealing because it estimates a distribution over functions, providing coherent estimates of uncertainty to map population variation. However GPR also has limitations: it is computationally prohibitive for large datasets and relies on predefined kernels with restricted functional form. Moreover, in the original implementation, brain measures were regressed independently (i.e., in a mass-univariate manner), which does not capitalize on the rich spatial structure of neuroimaging data. This last problem can be addressed by using multi-task GPR (MT-GPR) [Bonilla et al., 2008] to jointly predict multiple brain measurements. However, applying MT-GPR to neuroimaging data is very computationally demanding because of the need to invert large covariance matrices across both space and subjects. Recently, a combination of low-rank approximations and Kronecker algebra was proposed to scale MT-GPR to whole brain neuroimaging data [Kia and Marquand, 2018, Kia et al., 2018], which reduces the computational complexity with respect to the number of tasks by one order of magnitude. However, this comes with restrictive assumptions that the spatial structures of the signal and noise can be expressed by sets of orthogonal basis functions. Furthermore, its times complexity still remains cubic with the number of samples which is not appropriate for applications on large clinical cohorts.

Neural processes (NP) [Garnelo et al., 2018, Garnelo et al., 2018] are latent variable models that bring all the advantages of deep learning (e.g., representation learning and computationally efficient training and prediction) to the stochastic process framework and can address the problems described above. In the NP framework, a distribution over functions is modeled by learning an approximation to a stochastic process. Here, we present an application of NP to multivariate normative modeling of clinical neuroimaging data. This provides three advantages: i) like GPR, NP provides the necessary estimates of predictive uncertainty at test time; ii) similar to MT-GPR, it provides the possibility of learning structured variation; and iii) unlike alternatives, it is computationally scalable without restrictive assumptions on the orthogonality of lower dimensional representations of data. To this end, we make four contributions: i) in a tensor Gaussian predictive process (TGPP) framework [Kia et al., 2018], we formally define mixed-effect models of neuroimaging data [Friston et al., 1999] as stochastic processes; ii) we show how NP can be employed for mixed-effect modeling; iii) we use the resulting NP-based mixed-effect model to estimate a normative model of a clinical functional magnetic resonance imaging (fMRI) dataset; iv) we provide an example application of the proposed deep normative modeling for detecting psychiatric disorders in a novelty detection setting. Our experimental results show that the proposed method more accurately identifies ADHD patients from healthy individuals compared to the GP-based normative modeling.

§ METHODS

In this text, we use respectively calligraphic capital letters, $\mathcal{A}$, boldface capital letters, $\mathbf{A}$, and capital letters, $A$, to denote tensors, matrices, and scalars. We use $\times_1$ to denote $1st$-mode tensor product. We denote the vertical vector which results from collapsing the entries of a tensor $\mathcal{A}$ into a vector with $vec(\mathcal{A})$. Notation $\left | . \right |$ is accordingly used to represents the determinant of a matrix or the size of a set.

§.§ Mixed-Effect Modeling of MRIs in the TGPP Framework

Consider a neuroimaging study with $N$ subjects and let $\mathbf{X} \in \mathbb{R}^{N \times D}$ denote the design matrix of $D$ covariates of interest for $N$ subjects. Let $\mathcal{Y} \in \mathbb{R}^{N \times T_1 \times T_2 \times T_3}$  represent a $4$-order tensor of MRI data for corresponding $N$ subjects with respectively $T_1$, $T_2$, and $T_3$ voxels in $x$, $y$, and $z$ axes. In the normative modeling setting, we are interested in finding the function $f:\mathbf{X} \to \mathcal{Y}$. Adopting the tensor Gaussian predictive process (TGPP) [Kia et al., 2018] for structured multi-way mixed-effect modeling of MRI data, we have:

\begin{eqnarray} \label{eq:TGPP}
\mathcal{Y} = f(\mathbf{X}) = \mathbf{X} \times_1 \mathcal{A} + \mathcal{Z} + \mathcal{E} \quad ,
\end{eqnarray}

where $\mathcal{A} \in \mathbb{R}^{D \times T_1 \times T_2 \times T_3}$ represents the fixed-effect across subjects that contains regression coefficients estimated by solving the following linear equations:

\begin{eqnarray} \label{eq:fixed_effect}
\hat{\mathcal{Y}}[:,i,j,k] = \mathbf{X} \mathcal{A}[:,i,j,k], \quad & for \quad i=1,\dots, T_1;  \quad j=1,\dots, T_2; \quad k=1,\dots, T_3.
\end{eqnarray}

In eq:TGPP, $\mathcal{Z} \in \mathbb{R}^{N \times T_1 \times T_2 \times T_3}$ is the random-effect that characterizes the spatially structured joint variations from the fixed-effect across individuals in different dimensions of MRIs; and $\mathcal{E} \in \mathbb{R}^{N \times T_1 \times T_2 \times T_3}$ is heteroscedastic noise. Assuming a tensor-variate normal distribution for $\mathcal{Y}$ and a zero-mean tensor-variate normal distribution for $\mathcal{Z}+\mathcal{E}$, we have:

\begin{eqnarray} \label{eq:TVND}
p(\mathbf{X},\mathcal{Y}) = \mathcal{TN}(\hat{\mathcal{Y}}, \mathbf{S}) =
\frac{\exp(-\frac{1}{2}vec(\mathcal{Y}-\hat{\mathcal{Y}})^\top \mathbf{S}^{-1} vec(\mathcal{Y}-\hat{\mathcal{Y}}))}{\sqrt{(2\pi)^{NT} \left | \mathbf{S} \right |^{NT}}}  ,
\end{eqnarray}

where $\mathbf{S} \in \mathbb{R}^{NT \times NT}$ ($T = T_1 \times T_2 \times T_3$) is the covariance matrix of $\mathcal{Z} + \mathcal{E}$. Intuitively, the distribution of the mixed-effect in the joint hypercubic space of clinical covariates and neuroimaging measures can be described as a multi-dimensional Gaussian distribution with $vec(\hat{\mathcal{Y}})$ and $\mathbf{S}$ respectively serving as its mean and covariance.

§.§ Mixed-Effect Models of MRI Data as Stochastic Processes

The primary aim of this section is to formally define the structured mixed-effect model in eq:TGPP as stochastic process. This will provide the ingredients to employ NP for learning characteristics of the covariance matrix of the random-effect and noise in eq:TVND, i.e., $\mathbf{S}$.

Let $(\Omega, \Phi, \rho)$ represent a complete probability space (see Oksendal, 2003 or Appendix <ref> for definitions) where $\Omega$ is a set of clinical covariates and their corresponding neuroimaging measures pairs for $N$ subjects (i.e., $\left | \Omega \right | = N$) and $\Phi$ is a $\sigma$-algebra on $\Omega$ that contains all possible subsets of $\Omega$. Here, $\rho:\Phi \to [0,1]$ represents a probability measure that quantifies the probability of occurrence for any entry in $\Phi$. In this setting, each mixed-effect function $f_i$ estimated on the $i$th entry of $\Phi$ is a random variable, i.e., a $\Phi$-measurable function from $\Omega$ to a Borel set in $\mathbb{R}^{N \times T_1 \times T_2 \times T_3}$. Therefore, parametrizing $f_i$ on different subsets of $\Omega$; and considering the exchangeability and consistency properties of mixed-effect models [McCullagh, 2005, Nie and Yang, 2005], $\mathcal{Y}_i=f_i(\mathbf{X}_i) \mid_{i=1}^{\left | \Phi \right |}$ can be defined as stochastic processes [Garnelo et al., 2018]. As a corollary, for the $i$th entry in $\Phi$, $\phi_i = (\mathbf{X}_i,\mathcal{Y}_i) \subset \Omega$ with $\left | \phi_i \right | = N_i < N$, the joint distribution $p(\mathbf{X}_i,\mathcal{Y}_i)$ can be considered as a marginal for a higher-dimensional joint distribution in eq:TVND. We exploit this property to frame the problem of mixed-effect modeling in the neural processes framework [Garnelo et al., 2018]. To this end, given a particular realization of the mixed-effect stochastic process $f_i$, the joint distribution in eq:TVND can be rewritten as:

\begin{eqnarray} \label{eq:joint_dist}
p(\mathbf{X},\mathcal{Y}) = \sum_{i=1}^{\left | \Phi \right |} p(f_i) \mathcal{TN}(\mathcal{Y} \mid f_i, \mathbf{S}).
\end{eqnarray}

In an NP paradigm (see Appendix <ref> for background information on NP), we parametrize the integration over all $f_i(\mathbf{X})$ on a lower dimensional ($Q << T$) Gaussian distributed global latent variable $\mathbf{Z} \in \mathbb{R}^{N \times Q} \sim \mathcal{N}(\mu, \Sigma)$ where $f(\mathbf{X})=g(\mathbf{X},\mathbf{Z})$, resulting the following generative model:

\begin{eqnarray} \label{eq:generative_model}
p(\mathbf{Z}, \mathcal{Y} \mid \mathbf{X}) = p(\mathbf{Z}) \mathcal{TN}(\mathcal{Y} \mid g(\mathbf{X},\mathbf{Z}), \mathbf{S}) \quad ,
\end{eqnarray}

where $g(\mathbf{X},\mathbf{Z})$ is a deep neural network that learns the behavior of the mixed-effect model in an amortized variational inference regime [Kingma and Welling, 2013, Gershman and Goodman, 2014]. To this end, following the procedure proposed by Garnelo et al., 2018 the first challenge is to induce stochasticity, for which we need to define `context' and `target' points. While target points refer to full available information (e.g., all pixels in an image), the context points are intended to represent some partial information about the target function (e.g., a subset of pixels in an image). In this work, in order to adapt the NP for the mixed-effect modeling, we advance the concepts of context/target points [Garnelo et al., 2018] to context/target functions (see Section <ref> for discussion). The idea is to reduce the difference between the distribution of random context functions from the target function by minimizing their Kullback-Leibler (KL) divergence in the latent space. In our application in order to learn the distribution of the mixed-effect model in eq:TGPP, i.e., target function, we propose to use the estimated $\hat{\mathcal{Y}}_C \in \mathbb{R}^{N \times M \times T_1 \times T_2 \times T_3}$ (using eq:fixed_effect) on $M$ randomly drawn subsets of the training set as context functions. Then, using the actual corresponding neuroimaging training samples as target functions, the following evidence lower-bound should be optimized:

\begin{eqnarray} \label{eq:ELBO}
\log{p(\mathcal{Y} \mid \mathbf{X}, \hat{\mathcal{Y}}_C)} \geq \mathbb{E}_{q(\mathbf{Z} \mid \mathbf{X},\mathcal{Y})} \left [\log{p(\mathcal{Y} \mid \mathbf{Z}, \mathbf{X})}+\log{\frac{q(\mathbf{Z \mid \mathbf{X},\hat{\mathcal{Y}}_C})}{q(\mathbf{Z \mid \mathbf{X},\mathcal{Y}})}} \right] \quad ,
\end{eqnarray}

where $q(\mathbf{Z \mid \mathbf{X},\mathcal{Y}})$ is the variational posterior of the global latent variable that is parametrized on an encoder $h(\mathbf{X}, \hat{\mathcal{Y}}_C)$. In fact in this setting, each context function is a linear component of the target function that roughly approximates a stochastic process $f_i$. Having enough samples of context functions, large enough $M$, we expect the distribution of context functions to get rich enough to explain non-linear characteristics of the target function (i.e., the mixed-effect $f_i$). fig:demo shows a simplified illustration of this scenario in a 2D space where fitting enough linear models on subsets of noisy observations provides an estimation of the distribution of a non-linear target function. Furthermore, by minimizing the KL term in eq:ELBO, it is expected that the global latent variable $\mathbf{Z}$ will learn characteristics of the variance structure of the random-effect and noise terms (the diagonal elements of $\mathbf{S}$) from the difference between the context and target functions (recall that $\mathcal{Y} - \hat{\mathcal{Y}} = \mathcal{Z} + \mathcal{E}$).

A schematic illustration on approximating the distribution of a non-linear target function (red curve), e.g., a mixed-effect, from the distribution of linear context functions (blue lines), e.g., fixed-effects, which are fitted on $M$ random subsets of noisy observations (circles).

§.§ Deep Normative Modeling using Neural Processes

Using NP in the TGPP framework brings all the advantages of deep learning methods (e.g., representation learning from structured data and computational efficiency) for modeling the multi-way structured variation in neuroimaging data. It has been shown that modeling such structured variation provides the possibility of accurate unsupervised stratification of psychiatric patients in the normative modeling paradigm [Kia and Marquand, 2018, Kia et al., 2018]. To this end, here we introduce deep normative modeling, which utilizes an NP-based mixed-effect modeling and involves following three steps:

* Encoding phase: where an encoder $h(\mathbf{X}, \hat{\mathcal{Y}}_C)$ is learned to transfer the covariates, $\mathbf{X}$, and the estimated fixed-effects on $M$ randomly drawn samples from the training set, $\hat{\mathcal{Y}}_C$, to the parameters of the global latent variable $\mathbf{Z}$. Here, to preserve the 3D MRIs structure in the TGPP framework, we propose to use 3D-convolutional neural network (3D-CNN) layers to first transfer the $\hat{\mathcal{Y}}_C$ to a lower dimensional representation of neuroimages $\mathbf{R}_{\hat{\mathcal{Y}}} \in \mathbb{R}^{N \times T'}$. Note that using a CNN architecture in NP complicates fusing $\mathbf{X}$ with $\hat{\mathcal{Y}}_C$ in the encoder. When using fully-connected layers in the encoder (for example in Garnelo et al., 2018), this fusion is simply performed by concatenation. However, considering inherent structural differences between $\mathbf{X}$ and $\hat{\mathcal{Y}}_C$ the concatenation is impossible when using a CNN architecture. Therefore, this concatenation is performed in the latent output space $\mathbf{R}_{\hat{\mathcal{Y}}}$ (see Section <ref> for discussion on its advantages). Then, fully connected (FC) layers can be used to derive a latent representation in the joint space of clinical covariates ($\mathbf{X}$) and neuroimages, $\mathbf{R} \in \mathbb{R}^{N \times T''}$. It is worthwhile to emphasize that in this architecture, the aggregation across $M$ context functions is implicitly done by the 3D-CNN layers as they are considered as $M$ input channels to the CNN. Finally, two separate FC layers are used to transfer $\mathbf{R}$ to the means ($\mu_{\mathbf{Z}} \in \mathbb{R}^{N \times Q}$) and standard deviations ($\sigma_{\mathbf{Z}} \in \mathbb{R}^{N \times Q}$) of $\mathbf{Z}$.

* Decoding phase: where a decoder $g(\mathbf{X},\mathbf{Z})$ is learned to transfer back the joint covariates-latent space to the neuroimaging data $\mathcal{Y}$. Fully connected and 3D inverse CNN (3D-ICNN) layers can be accordingly used to reconstruct MRIs in the original space.

* Normative modeling: let $\mathcal{Y}^* \in \mathbb{R}^{N^* \times T_1 \times T_2 \times T_3}$ to represent the reconstructed neuroimaging data by the decoder $g(\mathbf{X}^*,\mathbf{Z})$ on $N^*$ test samples. Following Marquand et al., 2016 (see Appendix <ref> for details), we then compute statistical maps describing the deviation for each individual subject from the normative model, referred to as normative probability maps (NPMs), denoted by $\mathcal{N} \in \mathbb{R}^{N^* \times T_1 \times T_2 \times T_3}$ where $\mathcal{N} = (\mathcal{Y} - \mathcal{Y}^*)/\sqrt{\mathcal{S}}$. Here, $\mathcal{S}$ represents the sum of epistemic and aleatoric uncertainties, which respectively describe uncertainty about the true model  parameters and inherent variation in the data [Kendall and Gal, 2017]. To be able to calculate the epistemic uncertainty in our NP model, we keep the dropout layers active at test time [Gal and Ghahramani, 2016]. In the context of mixed-effect modeling of neuroimaging data (in eq:TGPP), the aleatoric uncertainty is the byproduct of two factors: i) the across-subject variability which is captured via the covariance of the random-effect $\mathcal{Z}$; and ii) noise in the data which is captured via covariance of $\mathcal{E}$. In the proposed NP framework, these two sources of uncertainties are learned from data and are summarized in the distribution of the global latent variable $\mathbf{Z}$. Therefore, given a test example of clinical covariates $\mathbf{x}^* \in \mathbf{X}^*$, we calculate the associated aleatoric uncertainty by sampling from the distribution of $\mathbf{Z}$.

An example NP architecture for mixed-effect modeling of MRIs.

§ EXPERIMENTAL MATERIALS AND SETUP

In our experiments, we use the response inhibition (i.e., `stop signal') task from the UCLA Consortium for Neuropsychiatric Phenomics dataset [Poldrack et al., 2016]. Specifically, we use the `Go' contrast volumes derived from the pipeline in Gorgolewski et al., 2017).[Available at <https://openfmri.org/dataset/ds000030/>.] The data consist of 119 healthy subjects; and 49, 39, and 48 individuals with schizophrenia (SCHZ), attention deficit hyperactivity disorder (ADHD), and bipolar disorder (BIPL), respectively. We cropped the volumes to the minimal bounding-box of $49 \times 61 \times 40$ voxels ($T_1=49,T_2=61,T_3=40,T=119560$). In order to accommodate the optimization scheme in eq:ELBO for fMRI data, the values of voxels are independently projected to the uniform $[0,1]$ interval using a robust quantile transformation. For clinical covariates, we use 11 factors of Barratt impulsiveness scores [Patton et al., 1995] ($D=11$) as impulsivity is a well-known feature for multiple psychiatric disorders and is implicated in response inhibition [Moeller et al., 2001].

We use three layers of 3D-CNNs followed by an FC layer to project $\hat{\mathcal{Y}}_C$ to $\mathbf{R}_{\hat{\mathcal{Y}}}$. In each CNN layer, we alternate a 3D-convolutional layer, a batch normalization layer [Ioffe and Szegedy, 2015], an average pooling layer, and a leaky ReLU activation function [Xu et al., 2015] (with negative slope of $0.01$). Then, two FC layers are used to transfer the merged $\mathbf{R}_{\hat{\mathcal{Y}}}$ and $\mathbf{X}$ to the middle joint representation $\mathbf{R}$. A similar reverse architecture is used for the decoder $g(\mathbf{X},\mathbf{Z})$ to transfer back the $\mathbf{Z}$ to $\mathcal{Y}$ space. fig:NP_architecture depicts a schematic of the employed NP architecture with detailed hyperparameter descriptions. Due to the small sample size and illustrative purpose of our experiments, we did not optimize the architecture and its hyperparameters (e.g., number of layers, number and the size of filters, number of neurons, etc.). The ADAM optimizer [Kingma and Ba, 2014] with decreasing learning rate (from $10^{-2}$ to $10^{-5}$) is used for optimization in 100 epochs.

We compare the normative models derived by NP and scalable multi-task Gaussian process tensor regression (sMT-GPTR) [Kia et al., 2018], in terms of their accuracy in detecting healthy subjects from patients.[The implementation for sMT-GPTR is available at <https://github.com/smkia/MTNorm>.] In the sMT-GPTR case, we set the number of basis functions across xyz dimensions of data $5,10$ and $3,5$ for the signal and noise, respectively (as they produced the best results in the original study). We evaluate normative modeling accuracy in a novelty detection scenario where we first train a model on a random subset of majority healthy subjects (75 healthy, 5 SCHZ, 5 ADHD, and 5 BIPL) and then calculate NPMs on a test set of remaining healthy subjects and patients. $\sim 16 \%$ of cases are included in the training set in order to seemingly simulate the average prevalence of general mental disorders in a cohort [Consortium, 2004]. We emphasize that the model has no access to the diagnostic labels during the training phase and thus our novelty detection approach is completely unsupervised. As in Marquand et al., 2016, we use extreme value statistics to provide a statistical model for the deviations (see Appendix <ref> for more details). Specifically, we use a block-maximum approach on the top 1% values in NPMs and fit these to a generalized extreme value distribution (GEVD) [Davison and Huser, 2015]. Then for a given test sample and given the shape parameter of GEVD, we compute the value of the cumulative distribution function of GEVD as the probability of that sample being an abnormal sample [Roberts, 2000]. Given these probabilities and actual labels, we evaluate the area under the receiver operating characteristic curve (AUC) to measure the performance of the model. All steps (random sampling, modeling, and evaluation) are repeated 10 times in order to estimate the fluctuations of models trained on different training sets. In all these experiments, ordinary least squares are used to estimate the fixed-effect (eq:fixed_effect) on bootstrapped subsets of the training set.[The scripts for experiments are available at <https://github.com/smkia/DNM>.]

Comparison between novelty detection performances of normative models derived by sMT-GPTR (with different number of bases for signal and noise) and NPs (with different $M$).

§ RESULTS

fig:novelty_detection compares the AUC of normative models derived by sMT-GPTR and NP. While sMT-GPTR shows slightly better performance in detecting SCHZ patients, NP provides substantially higher accuracy for ADHD cases. The methods perform similarly for BIPL. Considering the fact that these differences in performance are consistent across different model parameters and repetitions, it can be concluded that sMT-GPTR and NP are capturing different characteristics of the underlying biology of impulsivity. Furthermore, the above chance-level detection rates of NP models in ADHD and BIPL confirm a successful application of the proposed NP-based mixed-effect modeling in unsupervised diagnostic prediction. The significance of these results are even more pronounced considering the difficulty of the problem where a supervised support vector machine classifier provides only a chance-level performance in ADHD and BIPL cases (SCHZ$=0.67 \pm 0.07$, ADHD$=0.46 \pm 0.03$, BIPL$=0.47 \pm 0.06$).[See Kia et al., 2018 for training and evaluation configurations in the supervised setting.] Another important observation in NP models is the ascending trend of the detection performance as the number of samples from the fixed-effect ($M$) increases. This is compatible with the consistency property of mixed-effects as stochastic processes.

fig:analysis(a) depicts the average difference in NPMs of patient groups from the healthy population for NP(20) model (see Appendix <ref> for supplementary results). Different patterns of deviations from one diagnosis to another shed light on their different underlying biological causes. For example, the sign of deviations changes from SCHZ to ADHD patients in many regions. To further explore the link between these deviations and the level of impulsivity, we computed the coefficient of determination ($R^2$) between the average NPMs in 9 anatomical brain areas and the first principal component of covariates across different diagnostic groups (see fig:analysis(b)). The results show significantly (Bonferroni corrected F-test p-values) greater association between impulsivity and deviations in temporal lobes in ADHD and SCHZ patients compared to healthy individuals. This observation is compatible with previous research on the structural and functional engagement of temporal lobes in SCHZ and ADHD [Suddath et al., 1989, Kobel et al., 2010].

(a) The average difference between NPMs of healthy subjects and patients for NP(20). (b) $R^2$ between the impulsivity and deviation from the normative model across different anatomical brain areas ($** \quad p<0.01$ and $* \quad p<0.1$).

§ DISCUSSION

§.§ Toward Multivariate Normative Modeling on Large Clinical Cohorts

Including spatial information in probabilistic modeling and extending the mass-univariate normative modeling to its multivariate alternative is computationally very expensive. For $N$ samples and $T$ tasks, the time complexity of MT-GPR is cubic with respect to the number of samples and tasks, $\mathcal{O}(N^3 T^3)$. Many efforts have been devoted in order to reduce the time complexity of MT-GPR for large output spaces (i.e., large $T$) low-rank approximation and properties of Kronecker product [Alvarez and Lawrence, 2009, Stegle et al., 2011, Rakitsch et al., 2013, Kia and Marquand, 2018, Kia et al., 2018].  However, for very large sample-size datasets (i.e., for large $N$ and especially when $N>>T$), their time complexity still remains cubic with respect to $N$ that limits their applications in normative modeling on recently available large clinical cohorts [Sudlow et al., 2015] (with $N \approx 10^4-10^5$). One possible remedy for this problem is to approximate the posterior distribution of a probabilistic model with hidden variables in the stochastic variational inference framework [Hoffman et al., 2013]. Alvarez et al., 2010 made the first effort in employing variational inference in MT-GPR by introducing the variational inducing kernels that achieves a linear time complexity with respect to $N$. Our proposed NP-based normative modeling also employs the variational inference scheme, therefore, its computational complexity remains linear with respect to the number of samples in both training and inference phases [Garnelo et al., 2018]. Furthermore, since the spatial information is incorporated using a CNN architecture, there is no need to compute the inverse covariance matrix for the output space. These two properties make this method very suitable for multivariate normative modeling on large clinical cohorts of high-dimensional neuroimaging data.

§.§ From Context/Target Points to Context/Target Functions

In order to learn the joint distribution in eq:generative_model over random functions rather than a single function, the evidence lower-bound in eq:ELBO is optimized by minimizing the KL divergence between variational posteriors over context and target points. In the original NP framework [Garnelo et al., 2018], the target points are defined as whole points in the full dataset, while the context points are defined as a subset of target points that represent a partial knowledge about the full dataset. For example in the case of MNIST dataset, a random subset of pixels in an image can be used for context points. The random selection of pixels in the context points provides the desired stochasticity behavior in the NP framework. In this study, we advance the concepts of context/target points to target/context functions where the idea is to learn the distribution of a non-linear mixed-effect function, i.e., the target function, from a set of linear fixed-effect functions, i.e., context functions, estimated on a random subset of subjects. This alternation is key to learning the characteristics of the variance structure of random-effect and noise via the global latent variable and consequently it is crucial for normative modeling.

§.§ Preserving Spatial Structures via Convolutional Neural Processes

In this study, CNN-based architectures are proposed for encoding and decoding operations in NP. This results in two main advantages especially for the applications on neuroimaging data. First, it provides the possibility of preserving spatially-structured information in MRI data. Second, the parameter sharing gifted by CNN substantially reduces the computational costs in training and inference when dealing with very high-dimensional neuroimaging data.

§.§ Related Work

Rad et al., 2018 used a convolutional autoencoder for unimodal deep normative modeling of human movements recorded by wearable sensors. They used dropout technique in order to evaluate the parameter uncertainty of the model. Our proposed NP-based approach extends their effort in applying deep architectures to normative modeling from two perspectives: i) it provides the possibility of bimodal normative modeling. This is more appropriate for clinical usages where we are generally interested in interpreting the association between clinical covariates and biological measures [Marquand et al., 2016]; ii) using the fully probabilistic NP regime, we are also capable to evaluate aleatoric uncertainties resulting from individual differences and noise in addition to the epistemic parameter uncertainty.

§ CONCLUSIONS

In this paper, we proposed a principled approach for estimating spatially structured mixed-effects in neuroimaging data using neural processes. We demonstrated normative modeling as a possible target application for NP-based mixed-effect modeling. Even though the main focus in this study was on neuroimaging data, our contribution in framing the popular mixed-effect modeling as stochastic processes is quite general and opens the door for a wide range of NP applications in different research areas. Moreover, the presented application of NP for deep normative modeling of clinical neuroimaging data brings the advantages of deep neural networks in representation learning to the applications in precision psychiatry. Finally, the computational efficiency of NP in the training and evaluation phases (provided by its reliance on the variational inference) overcomes the lack of computational tractability of the GP-based normative modeling approaches especially when applied to large cohorts of high-dimensional neuroimaging data. For a possible future direction, we consider applying the proposed deep normative modeling approach to a large clinical neuroimaging cohort.

[Alvarez and Lawrence, 2009]
Mauricio Alvarez and Neil D Lawrence.
Sparse convolved Gaussian processes for multi-output regression.
In Advances in neural information processing systems, pages
57–64, 2009.

[Alvarez et al., 2010]
Mauricio A Alvarez, David Luengo, Michalis K Titsias, and Neil D Lawrence.
Efficient multioutput Gaussian processes through variational
inducing kernels.
In International Conference on Artificial Intelligence and
Statistics, pages 25–32, 2010.

[Bonilla et al., 2008]
Edwin V Bonilla, Kian M Chai, and Christopher Williams.
Multi-task Gaussian process prediction.
In Advances in neural information processing systems, pages
153–160, 2008.

[Bzdok and Meyer-Lindenberg, 2018]
Danilo Bzdok and Andreas Meyer-Lindenberg.
Machine learning for precision psychiatry: Opportunities and
Biological Psychiatry: Cognitive Neuroscience and
Neuroimaging, 30 (3):0 223 – 230, 2018.
ISSN 2451-9022.

[Consortium, 2004]
The WHO World Mental Health Survey Consortium.
Prevalence, Severity, and Unmet Need for Treatment of Mental
Disorders in the World Health Organization World Mental Health Surveys.
JAMA, 2910 (21):0 2581–2590, 06 2004.
ISSN 0098-7484.
URL <https://doi.org/10.1001/jama.291.21.2581>.

[Davison and Huser, 2015]
Anthony C Davison and Raphaël Huser.
Statistics of extremes.
Annual Review of Statistics and Its Application, 20
(1):0 203–235, 2015.
URL <https://doi.org/10.1146/annurev-statistics-010814-020133>.

[Foulkes and Blakemore, 2018]
Lucy Foulkes and Sarah-Jayne Blakemore.
Studying individual differences in human adolescent brain
Nature neuroscience, page 1, 2018.

[Friston et al., 1999]
Karl J Friston, Andrew P Holmes, CJ Price, C Büchel, and KJ Worsley.
Multisubject fMRI Studies and Conjunction Analyses.
NeuroImage, 100 (4):0 385 – 396, 1999.
ISSN 1053-8119.

[Gal and Ghahramani, 2016]
Yarin Gal and Zoubin Ghahramani.
Dropout as a bayesian approximation: Representing model uncertainty
in deep learning.
In Maria Florina Balcan and Kilian Q. Weinberger, editors,
Proceedings of The 33rd International Conference on Machine Learning,
volume 48 of Proceedings of Machine Learning Research, pages
1050–1059, New York, New York, USA, 20–22 Jun 2016. PMLR.
URL <http://proceedings.mlr.press/v48/gal16.html>.

[Garnelo et al., 2018]
Marta Garnelo, Dan Rosenbaum, Christopher Maddison, Tiago Ramalho, David
Saxton, Murray Shanahan, Yee Whye Teh, Danilo Rezende, and S. M. Ali Eslami.
Conditional neural processes.
In Jennifer Dy and Andreas Krause, editors, Proceedings of the
35th International Conference on Machine Learning, volume 80 of
Proceedings of Machine Learning Research, pages 1704–1713,
Stockholmsmässan, Stockholm Sweden, 10–15 Jul 2018a. PMLR.
URL <http://proceedings.mlr.press/v80/garnelo18a.html>.

[Garnelo et al., 2018]
Marta Garnelo, Jonathan Schwarz, Dan Rosenbaum, Fabio Viola, Danilo J Rezende,
SM Eslami, and Yee Whye Teh.
Neural processes.
arXiv preprint arXiv:1807.01622, 2018b.

[Gershman and Goodman, 2014]
Samuel Gershman and Noah Goodman.
Amortized inference in probabilistic reasoning.
In Proceedings of the annual meeting of the cognitive science
society, volume 36, 2014.

[Gorgolewski et al., 2017]
Krzysztof J Gorgolewski, Joke Durnez, and Russell A Poldrack.
Preprocessed consortium for neuropsychiatric phenomics dataset
[version 2; referees: 2 approved].
F1000Research, 60 (1262), 2017.

[Hoffman et al., 2013]
Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley.
Stochastic variational inference.
J. Mach. Learn. Res., 140 (1):0 1303–1347,
May 2013.
ISSN 1532-4435.
URL <http://dl.acm.org/citation.cfm?id=2502581.2502622>.

[Ioffe and Szegedy, 2015]
Sergey Ioffe and Christian Szegedy.
Batch normalization: Accelerating deep network training by reducing
internal covariate shift.
In Francis Bach and David Blei, editors, Proceedings of the
32nd International Conference on Machine Learning, volume 37 of
Proceedings of Machine Learning Research, pages 448–456, Lille,
France, 07–09 Jul 2015. PMLR.
URL <http://proceedings.mlr.press/v37/ioffe15.html>.

[Kapur et al., 2012]
Shitij Kapur, Anthony G Phillips, and Thomas R Insel.
Why has it taken so long for biological psychiatry to develop
clinical tests and what to do about it?
Molecular psychiatry, 170 (12):0 1174, 2012.

[Kendall and Gal, 2017]
Alex Kendall and Yarin Gal.
What uncertainties do we need in bayesian deep learning for computer
In Advances in Neural Information Processing Systems, pages
5580–5590, 2017.

[Kia and Marquand, 2018]
Seyed Mostafa Kia and Andre Marquand.
Normative modeling of neuroimaging data using scalable multi-task
gaussian processes.
In Alejandro F. Frangi, Julia A. Schnabel, Christos Davatzikos,
Carlos Alberola-López, and Gabor Fichtinger, editors, Medical Image
Computing and Computer Assisted Intervention – MICCAI 2018, pages 127–135,
Cham, 2018. Springer International Publishing.
ISBN 978-3-030-00931-1.

[Kia et al., 2018]
Seyed Mostafa Kia, Christian F. Backmann, and Andre F. Marquand.
Scalable multi-task gaussian process tensor regression for normative
modeling of structured variation in neuroimaging data.
arXiv preprint arXiv:1808.00036, 2018.

[Kingma and Ba, 2014]
Diederik P Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014.

[Kingma and Welling, 2013]
Diederik P Kingma and Max Welling.
Auto-encoding variational bayes.
arXiv preprint arXiv:1312.6114, 2013.

[Kobel et al., 2010]
Maja Kobel, Nina Bechtel, Karsten Specht, Markus Klarhöfer, Peter Weber, Klaus
Scheffler, Klaus Opwis, and Iris-Katharina Penner.
Structural and functional imaging approaches in attention
deficit/hyperactivity disorder: Does the temporal lobe play a key role?
Psychiatry Research: Neuroimaging, 1830 (3):0
230 – 236, 2010.
ISSN 0925-4927.

[Marquand et al., 2016]
Andre F Marquand, Iead Rezek, Jan Buitelaar, and Christian F Beckmann.
Understanding heterogeneity in clinical cohorts using normative
models: beyond case-control studies.
Biological psychiatry, 800 (7):0 552–561,

[McCullagh, 2005]
Peter McCullagh.
Exchangeability and regression models.
Oxford Statistical Science Series, 33:0 89, 2005.

[Mirnezami et al., 2012]
Reza Mirnezami, Jeremy Nicholson, and Ara Darzi.
Preparing for Precision Medicine.
New England Journal of Medicine, 3660 (6):0
489–491, 2012.
PMID: 22256780.

[Moeller et al., 2001]
F. Gerard Moeller, Ernest S. Barratt, Donald M. Dougherty, Joy M. Schmitz, and
Alan C. Swann.
Psychiatric aspects of impulsivity.
American Journal of Psychiatry, 1580 (11):0
1783–1793, 2001.
URL <https://doi.org/10.1176/appi.ajp.158.11.1783>.
PMID: 11691682.

[Nie and Yang, 2005]
Lei Nie and Min Yang.
Strong consistency of mle in nonlinear mixed-effects models with
large cluster size.
Sankhyā: The Indian Journal of Statistics, pages 736–763,

[Oksendal, 2003]
Bernt Oksendal.
Stochastic differential equations: an introduction with applications,

[Patton et al., 1995]
Jim H. Patton, Matthew S. Stanford, and Ernest S. Barratt.
Factor structure of the barratt impulsiveness scale.
Journal of Clinical Psychology, 510 (6):0
768–774, 1995.

[Poldrack et al., 2016]
Russell A Poldrack, Eliza Congdon, William Triplett, KJ Gorgolewski,
KH Karlsgodt, JA Mumford, FW Sabb, NB Freimer, ED London, TD Cannon, et al.
A phenome-wide examination of neural and cognitive function.
Scientific data, 3:0 160110, 2016.

[Rad et al., 2018]
Nastaran Mohammadian Rad, Twan van Laarhoven, Cesare Furlanello, and Elena
Novelty detection using deep normative modeling for imu-based
abnormal movement monitoring in parkinson’s disease and autism spectrum
Sensors, 180 (10), 2018.
ISSN 1424-8220.
URL <http://www.mdpi.com/1424-8220/18/10/3533>.

[Rakitsch et al., 2013]
Barbara Rakitsch, Christoph Lippert, Karsten Borgwardt, and Oliver Stegle.
It is all in the noise: Efficient multi-task Gaussian process
inference with structured residuals.
In Advances in neural information processing systems, pages
1466–1474, 2013.

[Roberts, 2000]
S.J. Roberts.
Extreme value statistics for novelty detection in biomedical data
IEE Proceedings - Science, Measurement and Technology,
147:0 363–367(4), November 2000.
ISSN 1350-2344.

[Stegle et al., 2011]
Oliver Stegle, Christoph Lippert, Joris M Mooij, Neil D Lawrence, and Karsten M
Efficient inference in matrix-variate gaussian models with iid
observation noise.
In Advances in neural information processing systems, pages
630–638, 2011.

[Suddath et al., 1989]
Richard L Suddath, Manuel F Casanova, Terry E Goldberg, David G Daniel, John R
Kelsoe Jr, and Daniel R Weinberger.
Temporal lobe pathology in schizophrenia: a quantitative magnetic
resonance imaging study.
American Journal of Psychiatry, 1460 (4):0
464–472, 1989.
URL <https://doi.org/10.1176/ajp.146.4.464>.
PMID: 2929746.

[Sudlow et al., 2015]
Cathie Sudlow, John Gallacher, Naomi Allen, Valerie Beral, Paul Burton, John
Danesh, Paul Downey, Paul Elliott, Jane Green, Martin Landray, Bette Liu,
Paul Matthews, Giok Ong, Jill Pell, Alan Silman, Alan Young, Tim Sprosen, Tim
Peakman, and Rory Collins.
Uk biobank: An open access resource for identifying the causes of a
wide range of complex diseases of middle and old age.
PLOS Medicine, 120 (3):0 1–10, 03 2015.
URL <https://doi.org/10.1371/journal.pmed.1001779>.

[Williams and Rasmussen, 1996]
Christopher KI Williams and Carl Edward Rasmussen.
Gaussian processes for regression.
In Advances in neural information processing systems, pages
514–520, 1996.

[Wolfers et al., 2015]
Thomas Wolfers, Jan K. Buitelaar, Christian F. Beckmann, Barbara Franke, and
Andre F. Marquand.
From estimating activation locality to predicting disorder: A review
of pattern recognition for neuroimaging-based psychiatric diagnostics.
Neuroscience and Biobehavioral Reviews, 57:0 328 –
349, 2015.
ISSN 0149-7634.

[Wolfers et al., 2018]
Thomas Wolfers, Nhat Trung Doan, Tobias Kaufmann, and et al.
Mapping the heterogeneous phenotype of schizophrenia and bipolar
disorder using normative models.
JAMA Psychiatry, 750 (11):0 1146–1155, 2018.
URL <+ http://dx.doi.org/10.1001/jamapsychiatry.2018.2467>.

[Wolfers et al., 2019]
Thomas Wolfers, Christian F. Beckmann, Martine Hoogman, Jan K. Buitelaar,
Barbara Franke, and Andre F. Marquand.
Individual differences v. the average patient: mapping the
heterogeneity in adhd using normative models.
Psychological Medicine, page 1–10, 2019.

[Xu et al., 2015]
Bing Xu, Naiyan Wang, Tianqi Chen, and Mu Li.
Empirical evaluation of rectified activations in convolutional
arXiv preprint arXiv:1505.00853, 2015.

[Zabihi et al., 2018]
Mariam Zabihi, Marianne Oldehinkel, Thomas Wolfers, Vincent Frouin, David
Goyard, Eva Loth, Tony Charman, Julian Tillmann, Tobias Banaschewski,
Guillaume Dumas, Rosemary Holt, Simon Baron-Cohen, Sarah Durston, Sven
Bölte, Declan Murphy, Christine Ecker, Jan K. Buitelaar, Christian F.
Beckmann, and Andre F. Marquand.
Dissecting the heterogeneous cortical anatomy of autism spectrum
disorder using normative models.
Biological Psychiatry: Cognitive Neuroscience and
Neuroimaging, 2018.
ISSN 2451-9022.

§ BACKGROUNDS

§.§ Neural Processes

A neural process (NP) [Garnelo et al., 2018] provides a computational tool to learn the distribution over a set of functions from distributions over a set of datasets $\Phi$. Assuming the $i$th dataset in $\Phi$ to contain a set of $N_i$ input-output pairs $(\mathbf{X}_i,\mathbf{Y}_i)$ where $\mathbf{X}_i \in \mathbb{R}^{N_i \times D}$ and $\mathbf{Y}_i \in \mathbb{R}^{N_i \times T}$ and we have $f_i: \mathbf{X}_i \to \mathbf{Y}_i$. For sake of simplicity, we refer to all $\mathbf{X}_i$ and $\mathbf{Y}_i$ in $\Phi$ as $\mathbf{X}$ and $\mathbf{Y}$, respectively. The goal of NP is to learn the distribution of $f_i$s from $(\mathbf{X},\mathbf{Y})$ pairs in $\Phi$ via learning the distribution of a global latent variable $\mathbf{Z}$ in the variational inference framework. For the generative model of an NP we have:

\begin{eqnarray} \label{eq:NP}
p(\mathbf{Z}, \mathbf{Y} \mid \mathbf{X}) = p(\mathbf{Z}) p(\mathbf{Y} \mid g(\mathbf{X},\mathbf{Z})) = \mathcal{N}(\mathbf{Y} \mid g(\mathbf{X},\mathbf{Z}), \mathbf{S}) \quad ,
\end{eqnarray}

where $g(\mathbf{X},\mathbf{Z})$ is the decoder function and parametrized by a neural network and $\mathbf{S}$ is the covariance matrix in the output space. Intuitively, the latent variable $\mathbf{Z}$ is intended to learn the statistical characteristics of the distribution of $f: \mathbf{X} \to \mathbf{Y}$. Then, the following approximation of variational posterior distribution is used in order to perform the approximate inference in NP:

\begin{eqnarray} \label{eq:NP_variational}
q(\mathbf{Z} \mid \mathbf{X}, \mathbf{Y}) = \mathcal{N}(m(\uplus h(\mathbf{X},\mathbf{Y})),s(\uplus h(\mathbf{X},\mathbf{Y}))) \quad ,
\end{eqnarray}

where, $h(\mathbf{X},\mathbf{Y})$ is the encoder function that is parametrized on neural network, $\uplus$ is the aggregator operator (for example, mean), and $m(.)$ and $s(.)$ are neural networks that map the aggregated values to the mean and standard deviation of $\mathbf{Z}$. Using the approximate variational posterior distribution in eq:NP_variational, the evidence lower bound (ELBO) on the log marginal likelihood is derived as follows:

\begin{eqnarray} \label{eq:NP_ELBO1}
\log{p(\mathbf{Y} \mid \mathbf{X})} \geq \mathbb{E}_{q(\mathbf{Z} \mid \mathbf{X},\mathbf{Y})} \left [\log{p(\mathbf{Y} \mid \mathbf{Z}, \mathbf{X})}+\log{\frac{p(\mathbf{Z})}{q(\mathbf{Z} \mid \mathbf{X},\mathbf{Y})}} \right].
\end{eqnarray}

In the NP framework, in order to learn such a distribution over random functions rather than a single function we need to create a context set of $M$ datasets $\Lambda \subset \Phi$ each of which containing input-output context pairs $(\mathbf{X}_\Lambda,\mathbf{Y}_\Lambda)$. These datasets are intended to represent some partial information about the target function $f:(\mathbf{X},\mathbf{Y})$. Thus, eq:NP_ELBO1 can be rewritten as:

\begin{eqnarray} \label{eq:NP_ELBO2}
\log{p(\mathbf{Y} \mid \mathbf{X}, \mathbf{X}_\Lambda,\mathbf{Y}_\Lambda)} \geq \mathbb{E}_{q(\mathbf{Z} \mid \mathbf{X},\mathbf{Y})} \left [\log{p(\mathbf{Y} \mid \mathbf{Z}, \mathbf{X})}+\log{\frac{p(\mathbf{Z} \mid \mathbf{X}_\Lambda,\mathbf{Y}_\Lambda)}{q(\mathbf{Z} \mid \mathbf{X},\mathbf{Y})}} \right] \quad ,
\end{eqnarray}

where the prior $p(Z)$ is replaced by the conditional prior $p(\mathbf{Z} \mid \mathbf{X}_\Lambda,\mathbf{Y}_\Lambda)$. Considering the intractability of this conditional prior we can approximate it by $q(\mathbf{Z} \mid \mathbf{X}_\Lambda,\mathbf{Y}_\Lambda)$. Therefore we optimize the following lower-bound in order to learn the distribution of $\mathbf{Z}$:

\begin{eqnarray} \label{eq:NP_ELBO3}
\log{p(\mathbf{Y} \mid \mathbf{X}, \mathbf{X}_\Lambda,\mathbf{Y}_\Lambda)} \geq \mathbb{E}_{q(\mathbf{Z} \mid \mathbf{X},\mathbf{Y})} \left [\log{p(\mathbf{Y} \mid \mathbf{Z}, \mathbf{X})}+\log{\frac{q(\mathbf{Z} \mid \mathbf{X}_\Lambda,\mathbf{Y}_\Lambda)}{q(\mathbf{Z} \mid \mathbf{X},\mathbf{Y})}} \right].
\end{eqnarray}

§.§ Normative Modeling

Normative modeling provides a framework for statistical inference on how the biological brain readouts of each individual subject deviate from the norm of a large population [Marquand et al., 2016]. Given $\mathbf{X} \in \mathbb{R}^{N \times D}$ and $\mathbf{Y} \in \mathbb{R}^{N \times T}$ respectively as matrix of $D$ clinical covariates and $T$ biological brain measures for $N$ subjects, normative modeling is performed in three steps:

* finding a mapping function $f: \mathbf{X} \to \mathbf{Y}$ from clinical covariates to brain readouts. While a wide range of linear and non-linear models can be used for this mapping. However, since computing the normative probability maps (see the next step) is strongly depends on estimating the prediction uncertainties, Bayesian regression approaches are the best candidates for normative modeling.

* calculating `normative probability maps' (NPMs), $\mathbf{Z} \in \mathbb{R}^{N \times T}$, as follows:

\begin{eqnarray} \label{eq:NPM}
\mathbf{Z} = \frac{\mathbf{Y} - \hat{\mathbf{Y}}}{\sqrt{\mathbf{S}}} \quad ,
\end{eqnarray}

where $\hat{\mathbf{Y}}$ and $\mathbf{S}$ are prediction mean and uncertainty, respectively. NPMs can be used to localize brain-related abnormalities at the single subject level [Wolfers et al., 2018, Zabihi et al., 2018, Wolfers et al., 2019]. To ensure accurate estimation of the NPMs it is important to model different sources of variation in data and model.

* computing subject-level summary statistics using a block-maximum approach by averaging top $1\%$ values in NPM of each subject. These summary statistics across subjects can be used as inputs to a novelty detection algorithm for diagnosis purposes [Kia and Marquand, 2018, Kia et al., 2018, Rad et al., 2018].

§.§ Novelty Detection using Generalized Extreme Value Distribution

According to Marquand et al., 2016, we can fit a generalized extreme value distribution (GEVD) on normative summary statistics across subjects in order to compute the abnormality index for each subject. This abnormality index can be defined as the probability of each sample being an abnormal sample by computing the cumulative distribution function of the fitted GEVD [Roberts, 2000]. For a random variable $a \in \mathbb{R}$, the cumulative distribution function of the GEVD is defined as below [Davison and Huser, 2015]:

\begin{eqnarray} \label{eq:GEVD}
F(a) = \left\{\begin{matrix}
\exp{(-[1+\xi (a-\mu)/\sigma]^{-1/\xi})}, &\xi \neq 0 \\
\exp{(-\exp{([-(a-\mu)/\sigma])})}, &\xi = 0
\end{matrix}\right. ,
\end{eqnarray}

$\mu \in \mathbb{R}$ and $\sigma>0$ are respectively the location and scale parameters and $\xi \in \mathbb{R}$ is the shape parameter. Depending on whether $\xi<0$, $\xi=0$, or $\xi>0$ the GEVD follows the special cases of the Weibull, Gumbel, Fr$\acute{e}$chet distributions, respectively.

§ SUPPLEMENTARY DEFINITIONS

Here are some complementary definitions from general probability theory to understand better the concepts in Section <ref>. The definitions are restated from Oksendal, 2003.

If $\Omega$ is a given set, then a $\sigma$-algebra $\Phi$ on $\Omega$ is a family $\Phi$ of subsets of $\Omega$ with the following properties:

* $\O \in \Phi$,

* $\forall \phi \in \Phi \Rightarrow \phi^C \in \Phi$, where $\phi^C$ is the complement set of $\phi$ in $\Omega$,

* $\phi_1, \phi_2, \dots \in \Phi \Rightarrow \bigcup_{i=1}^{\infty} \phi_i \in \Phi$.

Then, the pair $(\Omega,\Phi)$ is called a measurable space and the subsets of $\Omega$ that belong to $\Phi$ are called $\Phi$-measurable sets.

A probability measure $\rho$ on a measurable space $(\Omega,\Phi)$ is defined as a function $\rho:\Phi \to [0,1]$ such that:

* $\rho(\O) = 0, \rho(\Omega) = 1$,

* if $\phi_1, \phi_2, \dots \in \Phi$ and $\forall i, \forall j,i \neq j \Rightarrow \phi_i \cap \phi_j = \O$ then $\rho(\bigcup_{i=1}^{\infty} \phi_i)=\sum_{i=1}^{\infty} \rho(\phi_i)$.

The triple $(\Omega,\Phi,\rho)$ is called a probability space.

A probability space $(\Omega,\Phi,\rho)$ is called a complete probability space if $\Phi$ contains all subsets $\Lambda$ of $\Omega$ with P-outer measure zero, i.e., $\forall \phi \in \Phi$ with $\rho(\phi)=0$ we have $\forall \lambda \subset \phi \Rightarrow \lambda \in \Phi$. Any probability space can be made complete simply by adding to $\Phi$ all sets of outer measure 0 and by extending $\rho$ accordingly.

A stochastic process is a parametrized collection of random variables defined on a probability space $(\Omega,\Phi,\rho)$ and assuming values in $\mathbb{R}^n$.

§ SUPPLEMENTARY RESULTS

The average difference between NPMs of healthy subjects and patients for NP models with different $M$.

The average difference between NPMs of healthy subjects and patients for sMT-GPTR models with different number of basis functions for the signal and noise.