\section{Introduction}
\label{intro}
The time-dependent approach of an isolated quantum many-body system from an initial nonequilibrium state to thermal equilibrium
is an outstanding problem in many areas of physics. Examples at opposing ends of the energy scale are the local equilibration of quarks and gluons in the initial stages of relativistic heavy-ion collisions \cite{ba01}, and the thermalization of ultracold atoms, which can be accompanied by the time-dependent formation of a Bose--Einstein condensate (BEC) \cite{miesner_bosonic_1998}. Investigations of thermalization in relativistic heavy-ion collisions at energies reached at the Large Hadron Collider (LHC) have shown that the system has enough time to equilibrate locally. These results emerge from effective kinetic theories that rely on the quantum Boltzmann equation \cite{ba01,km11,jpb12,kur14,jpbl15,blmt17,fuku17}.

As a complementary mesoscopic approach, a nonlinear diffusion model has been developed \cite{gw18,bgw19,gw20,rgw20,sgw21,gw22} and is now used to calculate closed-form time-dependent fermionic and bosonic occupation-number distributions. The local equilibration of quarks and gluons in high-energy hadronic collisions, as well as the thermalization of bosonic ultracold atoms with time-dependent BEC formation is accounted for. A direct comparison with recent condensate-formation data for $^{39}$K from Cambridge University \cite{gli21} for various scattering lengths underlines the validity and usefulness of the analytical model for both infrared and ultraviolet dynamics.

The aim is to apply the nonlinear diffusion model \cite{gw18,gw22} in the high-energy regime to the thermalization of quarks and gluons as the relevant partonic constituents in relativistic heavy-ion collisions, and -- in the low-energy region -- to the thermalization of ultracold atoms, including the simultaneous time-dependent condensate formation. Only in the latter case is a direct comparison with experimental data possible, and it will be performed using results \cite{gli21} for ultracold potassium for various scattering lengths.
The kinetic nonlinear diffusion model will thus account for thermalization in both, fermionic and bosonic systems with the proper Fermi--Dirac and Bose--Einstein equilibrium limits. It is solved exactly through a nonlinear transformation, and applied to the local equilibration with suitably adapted transport coefficients in the high- and low-energy regions. The analytical model is complementary to, but simpler than existing numerical solutions of quantum Boltzmann equations.

Substantial efforts have been made to determine whether a transient gluon condensate can be formed in the initial stages of relativistic collision through elastic number-conserving gluon scatterings \cite{jpb12,jpb13}, bearing similarity to the BEC-formation in ultracold atoms below the critical temperature. It turns out, however, that inelastic number-changing collisions are much faster \cite{blmt17}, hindering condensate formation. In the present work, only inelastic gluon scatterings will be considered.

In contrast, thermalization in ultracold atoms occurs predominantly via number-conserving elastic collisions, triggering Bose-Einstein condensate formation \cite{an95,ket95,hul95,hul97}, and offering the possibility to measure the detailed time-dependence of condensate formation \cite{miesner_bosonic_1998,kdg02,gli21}. Some of these measurements have already been compared with numerical models for the time evolution of the growing condensate \cite{gz97,bzs00,kdg02}. The time-dependent condensate formation of ultracold $^{39}$K atoms that has recently been measured {using}
adjustable scattering lengths {\cite{gli21}} will be calculated and interpreted based on the analytical results of the nonlinear diffusion model.

\section{The model}
\label{model}
The nonlinear diffusion  model is formulated using the ergodic approximation \cite{snowo89,svi91,setk95,lrw96,jgz97}: The collision term for the single-particle
occupation-number distribution depends only on the energy, and on time. The time evolution of the distribution function
$n\equiv n_\alpha \equiv \langle n(\epsilon_\alpha,t)\rangle$  is governed by a nonlinear partial differential equation \cite{gw18,gw22}
 \begin{equation}
\partial_t {n}_{t}=-\partial_\epsilon\bigl[v\,n\,(1\pm n)+n\,\partial_\epsilon D\bigr]+\partial_{\epsilon\epsilon}\bigl[D\,n\bigr]\,,
 \label{boseq}
\end{equation}
where the $+$ sign gives rise to the bosonic, the $-$ sign to the fermionic solution.
In this nonlinear diffusion equation, the drift term $v\,(\epsilon,t)$ accounts for dissipative effects, whereas $D\,(\epsilon,t)$ mediates diffusion of particles in the energy space.
The many-body physics is contained in these
transport coefficients, which depend on energy, time, and the second moment of the interaction.

{The equation accounts for the time-dependent
thermalization in the cloud with particle number $N_\mathrm{th}(t)$, whereas condensate formation via elastic collisions will be considered through particle-number conservation at each timestep,
$N_\mathrm{c }(t) = N_\mathrm{i} – N_\mathrm{th}(t)$. The effect of interactions between particle cloud and condensate is included indirectly through the choice of the mesoscopic
transport coefficients $D, v.$

Alternatively, one could consider the condensate via an additional delta-function term in Eq.\,(1), as proposed for quantum Boltzmann equations \cite{jpbl15}. The modified diffusion equation would then conserve the particle number by itself. However, since for cold quantum gases one is interested in the time-dependent particle number in the condensate, the present approach without delta-function term seems more straightforward.}

The stationary solution $n_\infty(\epsilon)$ of the nonlinear diffusion equation is attained for $t\rightarrow\infty$ \cite{sgw21,gw22}, it equals the Bose--Einstein/Fermi--Dirac equilibrium distribution $n_\mathrm{eq}^\pm(\epsilon)$
\begin{equation}
n_\infty^\pm(\epsilon)=n_\mathrm{eq}^\pm(\epsilon)=\frac{1}{e^{(\epsilon-\mu)/T}\mp 1}\,,
 \label{Bose--Einstein}
\end{equation}
provided the ratio $v/D$ has no energy dependence, requiring the limiting value
$\lim_{t\rightarrow \infty}[-v\,(\epsilon,t)/D\,(\epsilon,t)] \equiv 1/T$.
The chemical potential is $\mu\leq0$ in a finite Bose system, and $\mu=\epsilon_\mathrm{f}>0$ for a Fermi system. It appears as an integration constant in the solution of
Eq.\,(\ref{boseq}) for $t\rightarrow\infty$.
For an energy-dependent diffusion coefficient, the term $n\,\partial D/\partial \epsilon$  in Eq.\,(\ref{boseq})
ensures to attain the correct equilibrium limit for $t\rightarrow \infty$ \cite{gw22}.

In the simplified case of constant transport coefficients, the nonlinear diffusion equation
preserves during the time evolution the essential features of Bose--Einstein and Fermi-Dirac
statistics that are contained in the bosonic or fermionic quantum Boltzmann equation.
{For constant transport coefficients, Eq.\,(1)} is one of the few nonlinear partial differential equations with a clear physical meaning that can be solved exactly through a nonlinear transformation \cite{gw18,gw22}. The resulting solution becomes
\begin{eqnarray}
        n(\epsilon,t) =\pm T\, \partial_\epsilon\ln{\mathcal{Z}(\epsilon,t)} \mp \frac{1}{2}= \pm \frac{T}{\mathcal{Z}}\, \partial_\epsilon \mathcal{Z} \mp
         \frac{1}{2}\,,
    \label{net}
    \end{eqnarray}
 where the time-dependent partition function ${\mathcal{Z}(\epsilon,t)}$ obeys a linear diffusion equation which has a Gaussian Green's function $G(\epsilon,x,t)$.

The partition function can be written as an integral over the Green's function and an exponential function $F(x)$

which depends on
 the initial occupation-number distribution $n_\mathrm{i}$
according to
 \begin{eqnarray}
    F(x) = \exp\Bigl[ -\frac{1}{2D}\bigl( v x\pm 2v \int_0^x n_{{i}}(y)\,{d}y \bigr) \Bigr]\,.
       \label{ini}
\end{eqnarray}

With the free Green's function $G\equiv G_\mathrm{free}$ of the linear diffusion equation, the physically correct solution is obtained for fermions \cite{bgw19}, but not for bosons due to the singularity in the infrared (IR): One has to consider the boundary conditions at the singularity $\epsilon = \mu \leq 0$ \cite{gw20}. They can be expressed as
 \(\lim_{\epsilon \downarrow \mu} n(\epsilon,t) = \infty\) \,$\forall$ \(t\). One obtains a vanishing partition function at the boundary \( \mathcal{Z} (\mu,t) = 0\), and the energy range is restricted to  $\epsilon \ge \mu$. This requires a new Green's function \cite{gw20,rgw20}
that equals zero at \(\epsilon = \mu\) $\forall \,t$. It can be written as
\begin{eqnarray}
    {G}_\mathrm{b} (\epsilon,x,t) = G_\mathrm{free}(\epsilon - \mu,x,t) - G_\mathrm{free}(\epsilon - \mu,-x,t),
    \label{eq:newGreens}
\end{eqnarray}
    and the partition function with this boundary condition becomes
     \begin{eqnarray}
    \mathcal{Z}_\mathrm{b} (\epsilon,t)= \int_{\mu}^{+\infty} G_\mathrm{b} (\epsilon,x,t)\,F(x)\,{d}x\,,
    \label{eq:partitionfunctionZ}
    \end{eqnarray}
which is equivalent to
$\int_{0}^{+\infty} G_\mathrm{b} (\epsilon,x,t)\,F(x+\mu)\,{d}x$:
The function $F$ remains unaltered with respect to Eq.\,(\ref{ini}), but its argument is shifted by the chemical potential.
With the free Green's function for fermions, and the bounded Green's function for bosons, we calculate the time-dependent partition
function, which can be done analytically for sufficiently simple initial conditions, and obtain the occupation-number distributions via the basic nonlinear transformation Eq.\,(\ref{net}).

\section{Thermalization of quarks}
 \begin{figure}
 \begin{center}
 \includegraphics[scale=0.6]{fig1.eps}%
\caption{\label{fig1} {{Thermalization of quarks.} Analytical solutions of the fermionic nonlinear diffusion equation for quarks
show the time-dependent approach to the Fermi--Dirac equilibrium distribution (solid curve) with $T=-D/v=600$\,MeV at times $t=1\times10^{-3}, 5\times10^{-3},  0.02, 0.05, 0.15~\mathrm{and}~0.4$ fm/$c$ (increasing dash lengths).
The transport coefficients are $D=1.44$ GeV$^2$/fm, $v=-2.4$ GeV/fm, the fermionic equilibration time is $\tau_\mathrm{eq}^\mathrm{F}=4D/v^2=1$ fm/$c$. The dotted curve is the Maxwell--Boltzmann distribution.}}
\end{center}
\end{figure}
The model is first applied to the thermalization of valence quarks in central relativistic heavy-ion collisions such as Pb--Pb at $\sqrt{s_\mathrm{NN}}=5.02$\,TeV. For an initial condition taken as a $\theta$-function up to the Fermi energy, $n_\mathrm{i}(\epsilon)=\theta({\epsilon_\mathrm{f}}-\epsilon)$, the function $F(x)$, the partition function $\mathcal{Z}$ and the single-particle distribution function can be evaluated analytically \cite{bgw19}. Typical results
for the transport coefficients $D, v$ of Table\,1 are shown in Fig.\,\ref{fig1} for six timesteps, a fermionic equilibration time $\tau^\mathrm{F}_\mathrm{eq}=4D/v^2=1$\,fm/$c$ as derived in \cite{gw18}, and an equilibrium temperature of $T=-D/v=600$\,MeV.

The time-dependent solutions at times beyond 0.1 fm/$c$ and the Fermi-Dirac equilibrium distribution at such a high initial temperature extend to negative energies, which is interpreted as antiparticle creation from the Fermi sea during the initial stages of the collision \cite{bgw19}. The initial temperature will then fall during expansion of the hot fireball until the hadronization temperature of about $T_\mathrm{H}\simeq160$\,MeV is reached, and light hadrons such as pions and kaons are formed.
\begin{table*}
\begin{center}
\hspace{4cm}
\caption{{Transport coefficients $D, v$ and equilibration times  $\tau_\mathrm{eq}$ for thermalization of quarks and gluons in central relativistic heavy-ion collisions, and of ultracold $^{39}$K atoms with Bose--Einstein condensate formation (scattering length $a=140\,a_0$).}}
\vspace{.2cm}
\label{tab1}
\begin{footnotesize}
\begin{tabular}{rrrrrr}
\hline\noalign{\smallskip}
&System~~~~~~&Energy~~~&$D$~~~~~~~~&$v$~~~~~~~~~~&$\tau_\mathrm{eq}$~~~~~\\
\noalign{\smallskip}\hline\noalign{\smallskip}
Quarks~~&$^{208}$Pb--$^{208}$Pb~~~~&5.02 TeV&1.44\,GeV$^2$/fm&$-2.4$\,GeV/fm~~&1.0\,fm/$c$\\
Gluons~~&$^{208}$Pb--$^{208}$Pb~~~~&5.02 TeV&1.20\,GeV$^2$/fm&$-2.0$\,GeV/fm~~&~~0.13\,fm/$c$\\
Atoms~~&$^{39}$K~~~~~~~~~~&195 nK~~&~~0.08\,nK$^2$/ms~~&$-0.00246$\,nK/ms&600 ms~~~\\
\noalign{\smallskip}\hline
\end{tabular}
\end{footnotesize}
\end{center}
\end{table*}
 \begin{figure}
 \begin{center}
 \includegraphics[scale=0.57]{fig2.eps}%
\caption{\label{fig2} {{Thermalization of gluons.} Analytical solutions of the nonlinear boson diffusion equation (NBDE) for gluons
with boundary conditions at the singularity show the time-dependent approach to the Bose--Einstein equilibrium distribution (solid curve) with $T=600$ MeV at times $t=2\times10^{-5},2\times10^{-4}, 2\times10^{-3}, 0.012, 0.04, 0.12~\mathrm{and}~0.4$ fm/$c$ (increasing dash lengths).
The transport coefficients are $D=1.2$ GeV$^2$/fm, $v=-2$ GeV/fm, the bosonic equilibration time is $\tau_\mathrm{eq}^\mathrm{B}=4D/(9v^2)=0.13$ fm/$c$. The dotted curve is the Maxwell--Boltzmann, the dot-dashed line the Rayleigh--Jeans distribution.}}
\end{center}
\end{figure}
\section{Thermalization of massless gluons}

To account for the local thermalization of gluons, the boundary conditions at the singularity must be considered. For massless gluons with
$\epsilon=p$ and $\theta$-function initial conditions
\cite{mue00} $n_\mathrm{i}=n_\mathrm{i}^0\,\theta(1-p/Q_\mathrm{s})$ with the gluon saturation momentum $Q_\mathrm{s}\simeq1$\,GeV and an average initial occupation $n_\mathrm{i}^0$, the function $F(x)$ can be evaluated in case of inelastic gluon scatterings with vanishing chemical potential $\mu=0$ as \cite{gw22}
\begin{eqnarray}
	\label{fxmu1}
		&F\,(x)=\exp[-v\,x/(2D)\\ \nonumber
		-&(v\,n_\mathrm{i}^0/D)(Q_\mathrm{s}-x)\,\theta\,(x-Q_\mathrm{s})+x]\,.

\end{eqnarray}
With the bounded Green's function Eq.\,(\ref{eq:newGreens}), the partition function and the gluon distribution function can be evaluated,
as shown in Fig.\,\ref{fig2} for $n_\mathrm{i}^0=1$ corresponding to an overoccupied situation \cite{jpb12} with the diffusion and drift coefficients of Table\,1. This results in a bosonic equilibration time $\tau^\mathrm{B}_\mathrm{eq}=4D/(9v^2)=0.13$\,fm/$c$ as derived in \cite{gw18} at $p=Q_\mathrm{s}$, and an equilibrium temperature of $T=-D/v=600$\,MeV.

The thermalization is faster as in case of quarks, because it is not hindered by Pauli's principle. The difference is enhanced by the colour factors $C_\mathrm{A}=3$ for gluons, and $C_\mathrm{f}=4/3$ for quarks.

{For gluons, the assumption of constant transport coefficients seems critical when viewed from the aspect of small-angle scattering equations \cite{jpb12,jpbl15,blmt17}.
However, when inserting physical timescales in the numerical calculations presented for inelastic gluon scattering in Ref.\,\cite{blmt17}, I obtain very similar analytical results for $D=\mathrm{const}$. There is presently no evidence that momentum-dependent diffusion coefficients are mandatory in gluon thermalization. }
 More detailed results including thermalization in an underoccupied system with the initial average occupation below the critical value \cite{jpb12} $r_\mathrm{c}=0.154$ for $\theta$-function initial conditions and inelastic gluon scatterings with $\mu=0$ are presented in Ref.\,\cite{gw22}.

\section{Thermalization of ultracold atoms}
 \begin{figure}
 \begin{center}
 \includegraphics[scale=0.58]{fig3.eps}%
\caption{\label{fig3} {{Thermalization of ultracold potassium atoms.}  The analytical solutions of the nonlinear boson diffusion equation (NBDE) for ultracold  $^{39}$K atoms are displayed as functions of time.
The initial state is a Bose--Einstein distribution with
a cut at $\epsilon_\mathrm{i}=15.56$\,nK, upper solid curve. The initial temperature is $T_\mathrm{i} = 130$\,nK, the final temperature $T_\mathrm{f} = -D/v = 32.5$\,nK (lower solid curve). Exact solutions of the NBDE are obtained as a finite series expansion with $k_\mathrm{max}=T_\mathrm{i}/T_\mathrm{f}=4$.
	The time evolution of the single-particle occupation-number distributions is shown at $t =8, 30, 100, 400$, and $4000$\, ms (increasing dash lengths). The transport coefficients are $D = 0.08$\,(nK)$^2$/ms,\,$v = -0.00246$\,nK/ms, the condensate formation time is $\tau_\mathrm{eq}^\mathrm{c}=f\,D/v^2=600$ ms, with a proportionality factor $f\simeq0.045$. The dot-dashed line is the Rayleigh--Jeans power law in the initial distribution.}}
	\end{center}
\end{figure}
\begin{figure}
\begin{center}
 \includegraphics[scale=0.57]{fig4.eps}%
\caption{\label{fig4} {{Time-dependent condensate fractions for an equilibrating 3D potassium Bose gas.}  The calculated condensate fractions based on the analytical solutions of the NBDE are displayed as functions of time (solid curves) and compared to Cambridge data \cite{gli21} (symbols, with error bars reflecting experimental fitting errors.). In the experiment, ultracold $^{39}$K atoms in a box trap are subjected to a rapid quench from $T_\mathrm{i} = 130$ nK to $T_\mathrm{f} = 32.5$ nK. Model results for an initial chemical potential
$\mu_\mathrm{i} = -0.67$ nK and a quench at $\epsilon_\mathrm{i} = 15.56$ nK are shown in a double-log plot for scattering lengths $a/a_0=$140 (circles), 280 (triangles), 400 (diamonds) and 800 (squares). The transport coefficients given in Table\,1 for
$a/a_0=140$ are scaled with the scattering lengths.}}
\end{center}
\end{figure}

Turning to the thermalization of ultracold atoms via number-conserving elastic collisions, the initial nonequilibrium distribution that is encountered in evaporative cooling --- or following a sudden energy quench \cite{gli21} --- can be represented as a Bose--Einstein equilibrium distribution with an energy cutoff at $\epsilon=\epsilon_\mathrm{i}$. The solution of the nonlinear boson diffusion equation (NBDE) for the combined initial- and boundary value problem \cite{gw20,rgw20} can still be calculated analytically, with the time-dependent partition function expressed as an infinite series
\begin{eqnarray}
{\mathcal{Z}}(\epsilon,t) &= \sqrt{4 D t} \, \exp\Bigl(-\frac{\mu}{2 T_\mathrm{f}}\Bigr)\qquad\qquad\\ \nonumber
 &\times\sum_{k=0}^{\infty} \dbinom{\frac{T_\mathrm{i}}{T_\mathrm{f}}}{k} (-1)^k \times
f_k^{{T}_\mathrm{i},{T}_\mathrm{f}} \,(\epsilon,t)\,.
\label{partfct}
\end{eqnarray}
The analytical expressions for $f_k^{{T}_\mathrm{i},{T}_\mathrm{f}} \,(\epsilon,t)$ are combinations of exponentials and error functions, they are given explicitly in
 Ref.\,\cite{rgw20}, and similarly, for the derivative of ${\mathcal{Z}}(\epsilon,t)$. The series terminates at $k_\mathrm{max}$ in case $T_\mathrm{i}/T_\mathrm{f}$ is an integer, which I use here to compute the exact solution for the thermalization and condensate formation of ultracold $^{39}$K atoms that has recently been investigated experimentally in Ref.\,\cite{gli21}, where $k_\mathrm{max}=T_\mathrm{i}/T_\mathrm{f}=4$.

The resulting time-dependent distribution functions following a quench at $\epsilon_\mathrm{i}=15.56$ nK of an initial thermal distribution of ultracold $^{39}$K atoms with $T_\mathrm{i}=130$\,nK corresponding to a removal of 77\% of the particles, 97.5\% of the total energy \cite{gli21}, and a final temperature $T_\mathrm{f}=32.5$\,nK that is below the critical temperature are shown in Fig.\,\ref{fig3} at five timesteps between 8\,ms and 4000\,ms, with transport coefficients $D, v$ from Table\,1 and an equilibration time of 600\,ms, which equals the condensate formation time: The condensate fraction reaches its equilibrium value within this time. The calculations are done for an interaction strength that corresponds to a scattering length $a=140\,a_0$, with $a_0=0.529\times10^{-10}$\,m the Bohr radius.
Details are presented in Ref.\,\cite{kgw22}.

With the analytical solutions of the nonlinear boson diffusion equation, we can now calculate the time-dependent condensate fraction from particle-number conservation as $N_\mathrm{c}(t)/N_\mathrm{i}=1-N_\mathrm{th}(t)/N_\mathrm{i}$, where $N_\mathrm{i}$ is the initial particle number just after the quench, and $N_\mathrm{th}(t)$ the integral of the time-dependent solution of the NBDE once the chemical potential has reached the value $\mu=0$ such that condensate formation starts,
\begin{equation}
N_\mathrm{th}(t)=\int_0^\infty g(\epsilon)\,n(\epsilon,t)\,{d}\epsilon\,.
\label{Nth}
\end{equation}
For a threedimensional box trap as in the Cambridge potassium experiment \cite{gli21}, the density of states is $g(\epsilon)=g_0\sqrt{\epsilon}$, and with the analytical nonequilibrium solutions $n(\epsilon,t)$ of Eqs.\,(\ref{net}),(\ref{partfct}), the integral is evaluated numerically
since presently no analytical solution is known, except for $t\rightarrow\infty$.

The Cambridge experiment has been performed for various scattering lengths, and here the analytical NBDE-results are compared in Fig.\,\ref{fig4} with data for $a/a_0=140, 280, 400$ and 800, where the transport coefficients start with the values given in Table\,1 for $a=140 a_0$, and are scaled with the scattering lengths rather than the cross sections due to the emerging coherence between the highly occupied IR-modes \cite{sto97,dwg17}.

{As in the high-energy case, the transport coefficients are assumed to be constant in the present cold-atom model calculations. An energy- or time-dependence of the diffusion coefficient would be required if the measured start of condensate formation was gradual, and not abrupt as is the case for constant diffusion coefficient.

There are indeed signs of a gradual onset in the historical $^{23}$Na MIT data \cite{miesner_bosonic_1998}, but not in the more recent $^{87}$Rb data, which have detailed error bars \cite{kdg02}. The current $^{39}$K deep-quench Cambridge data  \cite{gli21} show some indication of a gradual onset only at small scattering length $a=140a_0$  \cite{kgw22}.
Hence, the approximation $D=\mathrm{const}$ appears permissible also for cold atoms.}

For $a=140\,a_0$, data and model result almost coincide {beyond the onset}, whereas for the other scattering lengths, some deviations {larger than} the error bars are observed. In particular, the $1/a$-dependence of the timescales is not well fulfilled in the data for $a=400\,a_0$. But the overall agreement is rather satisfatory,
the equilibrium limit $N^\mathrm{eq}_\mathrm{c}/N_\mathrm{i}=1-g_0T_\mathrm{f}^{3/2}\zeta(3/2)\Gamma(3/2)/N_\mathrm{i}$ which equals the experimental value of $0.40(5)$ \cite{gli21} is approached within the condensate formation times.
The infrared power index $\alpha_\mathrm{IR}=1.08(9)$ that the experimentalists have extracted from their data is consistent with the NBDE calculation.
However, in the present formulation of the nonlinear diffusion model, thermalization proceeds continuously in time without nonthermal fixed points \cite{dwg17}.

These results confirm the validity and usefulness of the analytical nonlinear diffusion model and encourage further applications such as to the thermalization of fermionic cold atoms, but also in different areas of physics like reheating in inflationary cosmology.

\acknowledgments

Contributions of current or former Heidelberg students Thomas Bartsch, Johannes H\"olck, Anton Kabelac, Niklas Rasch and Alessandro Simon
to applications of the nonlinear diffusion model are gratefully acknowledged, as well as discussions with Thomas Gasenzer, Zoran Hadzibabic, Pascal Naidon, {and questions of the referee}.
I thank Emiko Hiyama for her hospitality at Tohoku University, Sendai, and RIKEN, Wako, where this work was completed with the support of {JSPS-BRIDGE fellowship BR200102}.

\input gw_epl_22.bbl

\end{document}


\section{Introduction}
\label{sec:intro}

These guidelines include complete descriptions of the fonts, spacing, and
related information for producing your proceedings manuscripts.

\section{Formatting your paper}
\label{sec:format}

All printed material, including text, illustrations, and charts, must be kept
within a print area of 7 inches (178 mm) wide by 9 inches (229 mm) high. Do
not write or print anything outside the print area. The top margin must be 1
inch (25 mm), except for the title page, and the left margin must be 0.75 inch
(19 mm).  All {\it text} must be in a two-column format. Columns are to be 3.39
inches (86 mm) wide, with a 0.24 inch (6 mm) space between them. Text must be
fully justified.

\section{Page title section}
\label{sec:pagestyle}

The paper title (on the first page) should begin 1.38 inches (35 mm) from the
top edge of the page, centered, completely capitalized, and in Times 14-point,
boldface type.  The authors' name(s) and affiliation(s) appear below the title
in capital and lower case letters.  Papers with multiple authors and
affiliations may require two or more lines for this information.

\section{Type-style and fonts}
\label{sec:typestyle}

To achieve the best rendering both in the printed and digital proceedings, we
strongly encourage you to use Times-Roman font.  In addition, this will give
the proceedings a more uniform look.  Use a font that is no smaller than nine
point type throughout the paper, including figure captions.

In nine point type font, capital letters are 2 mm high.  If you use the
smallest point size, there should be no more than 3.2 lines/cm (8 lines/inch)
vertically.  This is a minimum spacing; 2.75 lines/cm (7 lines/inch) will make
the paper much more readable.  Larger type sizes require correspondingly larger
vertical spacing.  Please do not double-space your paper.  True-Type 1 fonts
are preferred.

The first paragraph in each section should not be indented, but all the
following paragraphs within the section should be indented as these paragraphs
demonstrate.

\section{Major headings}
\label{sec:majhead}

Major headings, for example, "1. Introduction", should appear in all capital
letters, bold face if possible, centered in the column, with one blank line
before, and one blank line after. Use a period (".") after the heading number,
not a colon.

\subsection{Subheadings}
\label{ssec:subhead}

Subheadings should appear in lower case (initial word capitalized) in
boldface.  They should start at the left margin on a separate line.

\subsubsection{Sub-subheadings}
\label{sssec:subsubhead}

Sub-subheadings, as in this paragraph, are discouraged. However, if you
must use them, they should appear in lower case (initial word
capitalized) and start at the left margin on a separate line, with paragraph
text beginning on the following line.  They should be in italics.

\section{Printing your paper}
\label{sec:print}

Print your properly formatted text on high-quality, 8.5 x 11-inch white printer
paper. A4 paper is also acceptable, but please leave the extra 0.5 inch (12 mm)
empty at the BOTTOM of the page and follow the top and left margins as
specified.  If the last page of your paper is only partially filled, arrange
the columns so that they are evenly balanced if possible, rather than having
one long column.

In \LaTeX, to start a new column (but not a new page) and help balance the
last-page column lengths, you can use the command ``$\backslash$pagebreak'' as
demonstrated on this page (see the \LaTeX\ source below).

\section{Page numbering}
\label{sec:page}

Please do {\bf not} paginate your paper.  Page numbers, session numbers, and
conference identification will be inserted when the paper is included in the
proceedings.

\section{Illustrations, graphs, and photographs}
\label{sec:illust}

Illustrations must appear within the designated margins.  They may span the two
columns.  If possible, position illustrations at the top of columns, rather
than in the middle or at the bottom.  Caption and number every illustration.
All halftone illustrations must be clear black and white prints.  If you use
color, make sure that the color figures are clear when printed on a black-only
printer.

Since there are many ways, often incompatible, of including images (e.g., with
experimental results) in a \LaTeX\ document, below is an example of how to do
this \cite{Lamp86}.

\begin{figure}[htb]

\begin{minipage}[b]{1.0\linewidth}
  \centering
  \centerline{\includegraphics[width=8.5cm]{example-image}}
  \centerline{(a) Result 1}\medskip
\end{minipage}
\begin{minipage}[b]{.48\linewidth}
  \centering
  \centerline{\includegraphics[width=4.0cm]{example-image}}
  \centerline{(b) Results 3}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\includegraphics[width=4.0cm]{example-image}}
  \centerline{(c) Result 4}\medskip
\end{minipage}
\caption{Example of placing a figure with experimental results.}
\label{fig:res}
\end{figure}

\vfill
\pagebreak

\section{Footnotes}
\label{sec:foot}

Use footnotes sparingly (or not at all!) and place them at the bottom of the
column on the page on which they are referenced. Use Times 9-point type,
single-spaced. To help your readers, avoid using footnotes altogether and
include necessary peripheral observations in the text (within parentheses, if
you prefer, as in this sentence).

\section{Copyright forms}
\label{sec:copyright}

You must include your fully completed, signed IEEE copyright release form when
you submit your paper. We {\bf must} have this form before your paper can be
published in the proceedings.  The copyright form is available as a Word file,
a PDF file, and an HTML file. You can also use the form sent with your author
kit.

\section{Referencing}
\label{sec:ref}

List and number all bibliographical references at the end of the
paper.  The references can be numbered in alphabetic order or in order
of appearance in the document.  When referring to them in the text,
type the corresponding reference number in square brackets as shown at
the end of this sentence \cite{C2}.

\section{Compliance with ethical standards}
\label{sec:ethics}

IEEE-ISBI supports the standard requirements on the use of animal and
human subjects for scientific and biomedical research. For all IEEE
ISBI papers reporting data from studies involving human and/or
animal subjects, formal review and approval, or formal review and
waiver, by an appropriate institutional review board or ethics
committee is required and should be stated in the papers. For those
investigators whose Institutions do not have formal ethics review
committees, the principles  outlined in the Helsinki Declaration of
1975, as revised in 2000, should be followed.

Reporting on compliance with ethical standards is required
(irrespective of whether ethical approval was needed for the study) in
the paper. Authors are responsible for correctness of the statements
provided in the manuscript. Examples of appropriate statements
include:
\begin{itemize}
  \item ``This is a numerical simulation study for which no ethical
    approval was required.''
  \item ``This research study was conducted retrospectively using
    human subject data made available in open access by (Source
    information). Ethical approval was not required as confirmed by
    the license attached with the open access data.''
    \item ``This study was performed in line with the principles of
      the Declaration of Helsinki. Approval was granted by the Ethics
      Committee of University B (Date.../No. ...).''
\end{itemize}

\section{Acknowledgments}
\label{sec:acknowledgments}

IEEE-ISBI supports the disclosure of financial support for the project
as well as any financial and personal relationships of the author that
could create even the appearance of bias in the published work. The
authors must disclose any agency or individual that provided financial
support for the work as well as any personal or financial or
employment relationship between any author and the sources of
financial support for the work.

Other types of acknowledgements can also be listed in this section.

Reporting on real or potential conflicts of interests, or the absence
thereof, is required in the paper. Authors are responsible for
correctness of the statements provided in the manuscript. Examples of
appropriate statements include:
\begin{itemize}
  \item ``No funding was received for conducting this study. The
    authors have no relevant financial or non-financial interests to
    disclose.''
  \item ``This work was supported by […] (Grant numbers) and
    […]. Author X has served on advisory boards for Company Y.''
  \item ``Author X is partially funded by Y. Author Z is a Founder and
    Director for Company C.''
\end{itemize}

\bibliographystyle{IEEEbib}

\section{Introduction}
\label{sec:intro}

\subsection{Motivations}
Out-of-distribution (OOD) images are defined as being far from the training dataset distribution. Deep Learning (DL) models tend to behave inconsistently for this type of inputs, making OOD image detection crucial to avoid hidden model deficiencies \cite{yang2021generalized}. It is especially required in real-world automated pipelines, where input images may not be visually inspected before running the analysis. In the context of computerized medical-images analysis, a large variety of phenomenons in the input images can disturb a DL model and lead to unpredictable behaviours: noise, artifacts, variations in the imaging acquisition protocol and device, or pathological cases that were not included in the training dataset.

\subsection{Related Works}
OOD image detection frameworks can roughly be divided into two categories \cite{berger2021confidence}: methods that build a model specifically dedicated to the OOD task, and methods that rely on the prediction uncertainty of a trained model to detect abnormal inputs.

Within the first category, the most straightforward manner is to build a classifier directly on input images for the detection of OOD ones. For this, a Convolutional Neural Network (CNN) can be trained in a supervised manner, using an annotated dataset containing various types of real-world OOD \cite{bottani2022automatic}. Alternatively, Unsupervised Anomaly Detection (UAD) proposes to model the appearance of normal images by training an Auto-Encoder network (AE) to reconstruct in-distribution (ID) samples. At test-time, reconstruction is expected to be degraded for OOD samples, allowing for their detection \cite{gong2019memorizing,park2020learning}. Finally, in the context of semantic image segmentation, the Image Resynthesis framework \cite{lis2019detecting,xia2020synthesize} proposes to train a Generative Adversarial Network (GAN) to reconstruct the input image from the predicted segmentation. As the automated segmentation of OOD images is expected to be poor, the resynthesized input image should be distant from the true (unseen) input image.

While the former methods require the training of an auxiliary model (CNN, AE, or GAN) dedicated to the OOD detection task, uncertainty quantification (UQ) approaches propose to detect OOD inputs directly from the outputs of the trained predictive neural network. They rely on the hypothesis that the uncertainty of the deployed model should be higher when confronted to an unusual image, compared to the one of ID images. Based on this principle, OOD detection can be performed by monitoring the output uncertainty of the model. Various frameworks have been proposed to quantify uncertainty, and we present four of the most popular approaches below.

The baseline approach to estimate the uncertainty of a trained prediction model is to consider its output softmax probabilities. Following the intuition that the higher the probability, the more certain is the model, the \emph{Maximum Softmax Probability} (MSP) can be used as an estimator of uncertainty \cite{hendrycks2016baseline}. Another standard UQ methodology consists in producing a set of diverse and plausible predictions for the same input image. To achieve that, the \emph{MC dropout} approach \cite{gal2016dropout} consists in training a segmentation model with dropout and keep it activated during inference. At test time, each input image is forwarded multiple times through the model, with randomly sampled Dropout masks, yielding to multiple predictions for the same input image. Alternatively, \emph{Deep Ensemble} (DE) \cite{lakshminarayanan2017simple} proposes to train the same neural network multiple times. As the weights are randomly initialized at each training, the models reach different optimums. At test time, each individual model processes the input image independently, which produces a set of predictions. For both the MC dropout model and DE, uncertainty can then be estimated by computing the variance among the predictions. Finally, \emph{Deterministic Uncertainty Methods} (DUM) analyze the intermediate feature maps of a trained model to detect OOD inputs \cite{postels2021practicality}. It is based on the assumption that the hidden activations of the model should be different between an ID sample and an OOD one. In \cite{karimi2020improving}, authors propose to analyze the signature of an input image in the latent space of a segmentation model to detect OOD samples. At test time, an image OOD score is obtained by computing the distance between the signature of the test image, and the signatures of all training images.

As mentioned, OOD detection in medical image segmentation is critical, particularly in a clinical setting. Although uncertainty-based OOD detection has been applied for medical image segmentation \cite{karimi2020improving,mehrtash2020confidence,diao2022unified}, its evaluation is usually limited to extreme OOD cases (\textit{e.g.} different imaging modality, organ of interest or imaging center). While it is necessary to robustly detect such cases, we believe that the detection of more subtle artifacts (\textit{e.g.} image artifacts, truncation or skipped preprocessing) is also crucial in clinical practice.

\subsection{Contributions}
In this work, we propose to further explore uncertainty-based OOD detection in the context of medical image segmentation. As a use-case, we study OOD detection in input MRIs in the context of Multiple Sclerosis (MS) lesions segmentation. We develop a robust evaluation framework containing 14 synthetic or real abnormalities, representative of the abnormalities encountered in clinical routine. We use this framework to compare four of the most popular UQ-based OOD-detection models. Each model, initially binary, is additionally derived as a multiclass model that also segments anatomical regions of the brain in parallel with the MS-lesion segmentation. Indeed, we hypothesize that binary models trained to solely segment lesions tend to discard information not relevant for the task at hand, which may leave out key information for OOD detection. To verify this intuition, we compare the binary and multiclass versions of each segmentation model with respect to the OOD detection task.
We demonstrate that the retention of all additional information in multiclass models highly improves the OOD performance of uncertainty-based techniques, as compared to binary models.

\begin{figure*}[!t]
\centering
\includegraphics[width=\textwidth]{isbi_ood_grid.pdf}
\caption{In and Out-of-distribution data used in the experiments. In black: ID data and associated ground truth; blue: synthetic OOD; red: real OOD.}\label{fig:ood}
\end{figure*}

\section{Materials}

Data were collected to train the segmentation models, to evaluate the OOD detection, as well as to evaluate the segmentation performance of the models.

\subsection{In-distribution data}\label{ID_Data}
We collect data from several open-source MS datasets: the MSSEG-1 Challenge \cite{commowick2016msseg}, the WMH 2017 Challenge \cite{kuijf2019standardized}, and the database gathered by the University of Ljubljana, Slovenia (MSLUB) \cite{lesjak2018novel}. For each of these datasets, T1-weighted and T2-weighted FLAIR MRI sequences are provided, as well as the ground truth segmentation of White Matter Hyperintensities (WMH).
From the total of 143 subjects, we randomly set aside 50 subjects. The remaining 93 subjects are used to train the segmentation models while the 50 set-aside subjects are included in the OOD-detection evaluation database. These 50 subjects sharing the same origin than the 93 training images, they are considered as ID images.

We first co-register FLAIR and T1 sequences, followed by a rigid registration of both sequences on a 170 × 204 × 170 template with an isotropic resolution of 1mm. We use FastSurfer on the T1 sequences in order to obtain an anatomical atlas of the brain composed of 95 classes for each patient \cite{henschel2020fastsurfer}, which we further simplify into 6 brain regions: supratentorial white and grey matter, cerebrospinal fluid, infratentorial white and grey matter, and brain stem. We skull-strip the MRIs using the computed brain mask. We finally superpose the ground truth WMH delineations to the anatomical labels to obtain the final training masks composed of 7 foreground classes (6 anatomical classes plus 1 lesion class). Only the FLAIR sequences are used for the training of the segmentation models.

\subsection{Out-of-distribution data}\label{ood_data}
OOD data were collected or generated to build the OOD detection evaluation database (Figure \ref{fig:ood}). The same data were also used to evaluate the segmentation performance of the models, together with the 50 set-aside ID samples.

\subsubsection{Synthetic Out-of-distribution data}\label{syn_ood}
A common type of OOD in MRI consists is artifacts. However, finding real images with a controlled amount of artifacts to allow evaluation of OOD detection methods is difficult. We therefore generate realistic synthetic artifacted images from the set of ID evaluation images. We used the TorchIO Data Augmentation library \cite{perez2021torchio} to generate \emph{Downsample}, \emph{Bias}, \emph{Motion}, \emph{Spikes}, \emph{Gaussian Noise}, \emph{Ghost}, and \emph{Scaling} artifacts. We also implement a \emph{Truncation} operation that randomly removes the top and bottom slices of the MRI to mimic errors in the file transfer. From each ID evaluation image, we generate one OOD image per type of synthetic artefact. Having implemented 8 types of synthetic artifacts, we obtain 400 synthetically artifacted images.

\subsubsection{Real Out-of-distribution data}\label{real_ood}
We also propose several real-world OOD scenarios to validate our approaches. We use the set-aside T1-weighted sequences of the 50 test subjects, as well post-contrast T1-weighted (T1-Gd) sequences from the MSSEG-1 dataset to simulate mistakes in the input MRI sequences of an automated pipeline (T1 or T1-Gd sequences used as input, instead of expected FLAIR series). We use 50 FLAIR sequences from the BraTS 2020 dataset \cite{menze2014multimodal} containing brain tumors to test how the models behave on a pathology not observed during training. We also construct a version of the test dataset for which no skull-stripping is applied to mimic an error in the image preprocessing pipeline. Lastly, we use two other MS datasets that are not used during training to simulate domain-shift due to different acquisition protocols. We used 21 scans from the ISBI 2015 MS dataset \cite{carass2017longitudinal} and 60 scans from the MSSEG-2 dataset \cite{commowick2021msseg}. All data follow the same preprocessing steps presented in \ref{ID_Data}, except for the samples which are deliberately not skull-stripped.

\section{Method}
We train standard and MC dropout segmentation models in 2 settings: binary (with only WMH being segmented) and multiclass (WMH as well as the 6 brain anatomical regions are segmented). MSP is obtained from the softmax probabilities of the standard models, and DUM is implemented using the features of their penultimate layers \cite{karimi2020improving}. Each type of segmentation model is trained 5 times, allowing to obtain robust statistics for MSP, MC Dropout and DUM models, as well as to build DE models.

While DUM directly provides image-wise scores, the other implemented UQ methods (MSP, MC dropout and DE) provide one uncertainty score per voxel. For these methods, we compute the average of all the volume voxel uncertainties (estimated using the variance for MC dropout and DE) to get a single score per MRI \cite{karimi2020improving}.

We define OOD detection as a binary classification task, with ID images defined as negative samples, and OOD images as positive samples \cite{hendrycks2016baseline}. For each OOD dataset, we compute the image-level uncertainty scores and compare them with the scores obtained for the ID set of test images. We use the area under the ROC curve (AUROC) to estimate the ability to distinguish between ID and OOD images based on the image-level uncertainty scores. When the ground truth of WMH is available (i.e. not for the BraTS and MSSEG-2 data), we also report the WMH segmentation performances of MSP, MC dropout and DE, evaluated using the Dice score (DSC).

\begin{table*}[t!]
\resizebox{\textwidth}{!}{
\begin{tabular}{|l|cc|cc|cc|cc|cc|cc|c|c|}
\hline
& \multicolumn{4}{c|}{MSP} & \multicolumn{4}{c|}{MC dropout} & 	\multicolumn{4}{c|}{DE} & \multicolumn{2}{c|}{DUM} \\

 & \multicolumn{2}{c|}{Binary} & \multicolumn{2}{c|}{Multiclass} & \multicolumn{2}{c|}{Binary} & \multicolumn{2}{c|}{Multiclass} & \multicolumn{2}{c|}{Binary} & \multicolumn{2}{c|}{Multiclass} & Binary & Multiclass \\
  & AUROC  & DSC & AUROC  & DSC & AUROC & DSC & AUROC & DSC & AUROC & DSC & AUROC & DSC & AUROC & AUROC \\ \cline{2-15}

Test  & -  & 0.67 $\pm$ 0.00  & -   & 0.68 $\pm$ 0.00 & - & 0.66 $\pm$ 0.01 & - & 0.66 $\pm$ 0.01 & - & 0.68   & -  & 0.68 & -   & - \\

Downsample  & 0.50 $\pm$ 0.02  & 0.62 $\pm$ 0.00  & \textbf{0.94 $\pm$ 0.01} & 0.63 $\pm$ 0.00 & 0.52 $\pm$ 0.02 & 0.61 $\pm$ 0.01 & 0.84 $\pm$ 0.03  & 0.62 $\pm$ 0.01 & 0.53 & 0.63 & 0.88  & 0.64 & 0.88 $\pm$ 0.03 & \textbf{0.94 $\pm$ 0.02} \\

Bias  & 0.58 $\pm$ 0.11  & 0.55 $\pm$ 0.01  & 0.74 $\pm$ 0.02 & 0.59 $\pm$ 0.01 & 0.51 $\pm$ 0.04  & 0.50 $\pm$ 0.02 & 0.71 $\pm$ 0.02  & 0.58 $\pm$ 0.02 & 0.58 & 0.56 & 0.74 & 0.59 & \textbf{0.94 $\pm$ 0.01} & 0.92 $\pm$ 0.02 \\

Motion  & 0.29 $\pm$ 0.16  & 0.63 $\pm$ 0.00 & \textbf{0.99 $\pm$ 0.01} & 0.64 $\pm$ 0.01 & 0.40 $\pm$ 0.01  & 0.63 $\pm$ 0.01 & 0.91 $\pm$ 0.03 & 0.63 $\pm$ 0.01 & 0.46 & 0.64 & 0.96 & 0.65 & 0.98 $\pm$ 0.02 & \textbf{0.99 $\pm$ 0.00} \\

Spikes & 0.15 $\pm$ 0.16 & 0.19 $\pm$ 0.04  & \textbf{1.00 $\pm$ 0.00} & 0.15 $\pm$ 0.04 & 0.51 $\pm$ 0.14  & 0.31 $\pm$ 0.05 & \textbf{1.00 $\pm$ 0.00}  & 0.26 $\pm$ 0.07 & 0.57 & 0.18   & \textbf{1.00} & 0.18 & \textbf{1.00 $\pm$ 0.00} & \textbf{1.00 $\pm$ 0.00} \\

Noise & 0.40 $\pm$ 0.36  & 0.64 $\pm$ 0.00  & 0.98 $\pm$ 0.01    & 0.64 $\pm$ 0.01  & 0.38 $\pm$ 0.03  & 0.64 $\pm$ 0.01 & 0.93 $\pm$ 0.05  & 0.64 $\pm$ 0.02 & 0.57 & 0.65 & 0.97 & 0.65  & \textbf{1.00 $\pm$ 0.00} & 0.99 $\pm$ 0.00 \\

Ghost & 0.41 $\pm$ 0.20 & 0.65 $\pm$ 0.00 & 0.97 $\pm$ 0.03 & 0.65 $\pm$ 0.01 & 0.48 $\pm$ 0.02  & 0.64 $\pm$ 0.01 & 0.94 $\pm$ 0.03  & 0.64 $\pm$ 0.01 & 0.52  & 0.66 & 0.96 & 0.66 & 0.96 $\pm$ 0.01 & \textbf{0.99 $\pm$ 0.04} \\

Truncation & 0.56 $\pm$ 0.03  & 0.64 $\pm$ 0.00  & 0.15 $\pm$ 0.01  & 0.63 $\pm$ 0.00  & 0.52 $\pm$ 0.01  & 0.62 $\pm$ 0.01 & 0.49 $\pm$ 0.04  & 0.62 $\pm$ 0.01 & 0.52 & 0.65 & 0.67 & 0.65 & 0.89 $\pm$ 0.01 & \textbf{0.96 $\pm$ 0.01} \\

Scale & 0.48 $\pm$ 0.11  & 0.61 $\pm$ 0.00  & 0.95 $\pm$ 0.03 & 0.61 $\pm$ 0.01 & 0.59 $\pm$ 0.03  & 0.59 $\pm$ 0.01 & 0.95 $\pm$ 0.01  & 0.60 $\pm$ 0.01 & 0.58    & 0.62   & \textbf{0.99} & 0.63 & \textbf{0.99 $\pm$ 0.00} & \textbf{0.99 $\pm$ 0.00} \\

T1 & 0.46 $\pm$ 0.21  & 0.00 $\pm$ 0.00  & \textbf{1.00} $\pm$ 0.00  & 0.00 $\pm$ 0.00    & 0.78 $\pm$ 0.12  & 0.00 $\pm$ 0.00 & 0.99 $\pm$ 0.00  & 0.00 $\pm$ 0.00 & 0.86    & 0.00   & \textbf{1.00} & 0.00 & 0.96 $\pm$ 0.0  & 0.99 $\pm$ 0.00            \\

T1-Gd & 0.31 $\pm$ 0.20  & 0.00 $\pm$ 0.00  & \textbf{1.00} $\pm$ 0.00 & 0.00 $\pm$ 0.00 & 0.42 $\pm$ 0.12 & 0.00 $\pm$ 0.00 & 1.00 $\pm$ 0.00  & 0.00 $\pm$ 0.00 & 0.49 & 0.00 & \textbf{1.00} & 0.00 & 0.98 $\pm$ 0.01  & \textbf{1.00 $\pm$ 0.00} \\

BraTS & 0.54 $\pm$ 0.17  & - & \textbf{1.00} $\pm$ 0.00 & - & 0.93 $\pm$ 0.02  & - & \textbf{1.00 $\pm$ 0.00}  & - & 0.96   & - & \textbf{1.00} & - & 0.99 $\pm$ 0.00 & \textbf{1.00 $\pm$ 0.00} \\

Skull  & 0.35 $\pm$ 0.17 & 0.51 $\pm$ 0.01  & \textbf{1.00 $\pm$ 0.00} & 0.52 $\pm$ 0.02 & 0.58 $\pm$ 0.06 & 0.49 $\pm$ 0.01 & \textbf{1.00 $\pm$ 0.00} & 0.54 $\pm$ 0.02 & 0.68 & 0.53 & \textbf{1.00} & 0.56 & 0.97 $\pm$ 0.01 & 0.99 $\pm$ 0.00 \\

ISBI 2015 & 0.34 $\pm$ 0.20 & 0.62 $\pm$ 0.01 & \textbf{1.00 $\pm$ 0.00} & 0.65 $\pm$ 0.02 & 0.46 $\pm$ 0.05 & 0.65 $\pm$ 0.02 & 0.99 $\pm$ 0.01 & 0.65 $\pm$ 0.02 & 0.45 & 0.62 & \textbf{1.00}  & 0.67 & 0.97 $\pm$ 0.02 & 0.99 $\pm$ 0.00  \\

MSSEG-2 & 0.62$\pm$ 0.05  & - & 0.94 $\pm$ 0.01 & - & 0.51 $\pm$ 0.02 & - & 0.93 $\pm$ 0.04 & - & 0.50 & - & \textbf{0.99} & - & 0.95 $\pm$ 0.02 & \textbf{0.99 $\pm$ 0.00} \\ \hline

\end{tabular}}
\caption{AUROC and DSC scores obtained for each ID and OOD datasets. For all methods except DE, we present the performance averaged over the 5 runs, as well as the standard deviation. The best OOD detection performance for each dataset is highlighted in \textbf{bold}. }\label{tab:results}
\end{table*}

\section{Implementation details}
Our framework is implemented using PyTorch \cite{paszke2019pytorch}. We use a 3D Residual UNet \cite{kerfoot2018left} as the backbone of the segmentation models. For MC dropout, we apply a 3D dropout with a rate of $20\%$ in each layer of the encoder and decoder. The models are trained using 3D patches of $160\times192\times32$ with the Dice loss \cite{milletari2016v}. We use a simple Data Augmentation pipeline during training composed of flipping, rotation, and contrast alteration. At inference, we perform $T=20$ inferences per image to obtain the final segmentation and uncertainty estimates for MC dropout models.

\section{Results}
The AUROC and DSC scores for each OOD dataset and UQ method are presented in Table \ref{tab:results}. As expected, DSC scores are inferior on OOD data (\emph{e.g.} \emph{T1-Gd} dataset), meaning that the inability to detect OOD could result in catastrophic outcomes in the context of a real world automated pipeline. This is the case with binary implementations of MSP, MC dropout and DE that reach low OOD detection performance on most datasets, being even worst than chance ($\text{AUROC} < 0.5$) in some settings (\emph{e.g.}, MSP Binary on the \emph{Spikes} dataset). Exceptions are made for the MC dropout Binary and DE Binary that show satisfying OOD detection on the T1 and BraTS datasets, which are however the most obvious types of OOD. On the contrary, the multiclass versions demonstrate high improvement with respect to the OOD task, being able to perfectly distinguish between ID and OOD samples in many settings (\emph{e.g.}, $\text{AUROC}=1.00$ for DE Multiclass on the \emph{Skull} dataset). Interestingly, DUM models reach high detection accuracy in both settings, binary and multiclass.

\section{Discussion and Conclusion}
In this work, we investigate the use of UQ approaches to detect OOD inputs in the context of medical image segmentation. We elaborate an evaluation framework with a set of 14 various synthetic or real OODs to evaluate several popular UQ methods for OOD detection. Our experimental results demonstrate that standard UQ approaches such as MSP, MC dropout and DE fail to detect OOD inputs when applied on binary segmentation models, even on extreme OOD cases (\emph{e.g.} \emph{T1-Gd} dataset). Overall, binary segmentation models tend to be overconfident on abnormal inputs.

Second, our experiments show that UQ methods applied on models that were trained to segment anatomical labels additionally to lesions (multiclass) achieve excellent OOD detection. Intuitively, this additional segmentation task enforces the model to take into account the entirety of the available brain information to make its prediction, increasing the quantity of information available to decide whether the input is abnormal or not.

Finally, DUM models show excellent OOD detection performances, with DUM Multiclass being a little more consistent than the Binary version. This suggests that the hidden activations of a trained model are a powerful indicator of the distance of a sample to the training set. Interestingly, DUM is very computationally efficient, as they only require a single trained model, in contrast to the MC dropout and DE that require either performing multiple predictions with the same network, or training several models. These approaches are thus particularly interesting for the implementation of OOD-detection tasks in a clinical setting.

While multiclass seems beneficial for OOD detection, it requires the generation of anatomical atlases to train the model, which may not always be feasible. In that case, the DUM approach is the most suitable, as it achieves strong OOD detection results even on binary models.

Future work will investigate the impact of the number of learned anatomical labels on the OOD detection and compare the performance of uncertainty-based method to the performance of dedicated models such as AE or GANs. We will also investigate the impact of preprocessing choices (\emph{e.g.} skullstripping, debiasing) on OOD detection performance.

\section{Acknowledgments}
BL, AT, SD are employees of the Pixyl company. MD and FF serve on Pixyl scientific advisory board.

\section{Compliance with Ethical Standards}
This research study was conducted retrospectively using human subject data made available in open access. Ethical approval was not required as confirmed by the license attached with the open access data.

\bibliographystyle{IEEEbib}