# Breather gas fission from elliptic potentials in self-focusing media

Gino Biondini1,2, Gennady A. El3 Xu-Dan Luo4, Jeffrey Oregero5, and Alexander
Tovbis6
1 Department of Mathematics, State University of New York at Buffalo, Buffalo,
New York 14260, USA  2 Department of Physics, State University of New York at
Buffalo, Buffalo, New York 14260, USA  3 Department of Mathematics, Physics
and Electrical Engineering, Northumbria University, Newcastle upon Tyne, NE1
8ST, United Kingdom  4 Academy of Mathematics and Systems Science, Chinese
Academy of Sciences, Beijing 100190, China  5 Department of Mathematics,
University of Kansas, Lawrence, KS 66045, USA  6 Department of Mathematics,
University of Central Florida, Orlando, Florida 32816, USA

###### Abstract

We present an analytical model of integrable turbulence in the focusing
nonlinear Schrödinger (fNLS) equation, generated by a one-parameter family of
finite-band elliptic potentials in the semiclassical limit. We show that the
spectrum of these potentials exhibits a thermodynamic band/gap scaling
compatible with that of soliton and breather gases depending on the value of
the elliptic parameter $m$ of the potential. We then demonstrate that, upon
augmenting the potential by a small random noise (which is inevitably present
in real physical systems), the solution of the fNLS equation evolves into a
fully randomized, spatially homogeneous breather gas, a phenomenon we call
breather gas fission. We show that the statistical properties of the breather
gas at large times are determined by the spectral density of states generated
by the unperturbed initial potential. We analytically compute the kurtosis of
the breather gas as a function of the elliptic parameter $m$, and we show that
it is greater than 2 for all non-zero $m$, implying non-Gaussian statistics.
Finally, we verify the theoretical predictions by comparison with direct
numerical simulations of the fNLS equation. These results establish a link
between semiclassical limits of integrable systems and the statistical
characterization of their soliton and breather gases.

Introduction. The focusing nonlinear Schrödinger (fNLS) equation is a
ubiquitous model describing nonlinear wave propagation arising in a variety of
physical settings, including deep water waves [1], optics [2], plasmas [3],
and Bose-Einstein condensates (BECs) [4]. The fNLS equation is also a
completely integrable infinite-dimensional Hamiltonian system [5, 6, 7, 8, 9,
10], endowed with a deep mathematical structure, including the existence of
infinite families of exact solutions with both zero and non-zero background
describing the elastic interactions of $N$ solitons. In addition, its initial
value problem can in some cases be solved by the inverse scattering transform
(IST) [11, 12, 13, 14, 15, 16, 17, 18, 19]. The formulation of the IST is
based on the representation of the nonlinear evolution equation as the
compatibility condition of two linear equations called a Lax pair. The first
half of the Lax pair of the fNLS equation, namely the Zakharov-Shabat (ZS)
scattering problem, is equivalent to an eigenvalue problem for a non-self-
adjoint one-dimensional Dirac operator.

While the classical theory and applications of the fNLS equation are mostly
concerned with the description of regular, deterministic wave structures, the
inherent statistical nature of some physical wave phenomena in focusing media
(e.g., rogue wave emergence) calls for the study of stochastic fNLS solutions,
characterized in terms of the probability density function, correlation
function etc. Establishing a connection between the IST spectra of random fNLS
solutions and their statistical properties in physical space represents a
challenging problem, which has recently been formulated in the context of
integrable turbulence [20, 21] — the general theoretical framework for the
description of a broad spectrum of stochastic wave phenomena in physical
systems modeled by integrable equations. A particular type of integrable
turbulence, termed soliton gas (SG), has recently attracted considerable
attention [23, 24] due to its appearance in many physical systems including
water waves [25, 26, 27, 28], nonlinear optics [30, 29, 31] and BECs [32].
Notably, SG dynamics has been shown to underpin some fundamental physical
phenomena such as spontaneous modulational instability [33] and the rogue wave
emergence [34] in focusing media.

The concept of a SG was introduced in [35] as an infinite collection of
randomly distributed solitons with small spatial density and with a certain
amplitude distribution. Soliton interactions, accompanied by well-defined
phase shifts, result in a modification of the effective velocity of a “tracer”
soliton in a gas over large propagation distances, enabling an approximate
description of the emergent, large-scale hydrodynamics or kinetics of a weakly
non-uniform/non-equilibrium SG. The kinetic description of a SG was
generalized in [22] to a dense SG using finite-gap theory, and a general
phenomenological construction of SG kinetic equation for a broad class of
integrable systems was proposed in [36]. A systematic spectral theory of fNLS
SGs was developed in [37], where it was also extended to the case of solutions
on a non-zero background, i.e., a breather gas (BG). For recent advances on
the spectral theory of soliton and breather gases and its relation to the
generalised hydrodynamics of integrable many-body classical and quantum
systems, see [24].

While BGs represent a natural generalization of SGs, possible mechanisms for
their generation have remained largely unexplored. In this Letter we present
an analytically tractable model of BG generation via fission, based on the
semiclassical limit of the fNLS equation with initial data in the form of a
periodic elliptic “$\mathop{\rm dn}\nolimits$” potential with elliptic
parameter $m\in[0,1]$, augmented by a small random noise. It was shown in [38]
that the ZS spectrum of such elliptic potentials can be characterized
analytically. Here we show that, in the semiclassical limits, these potentials
are compatible with the so-called thermodynamic spectral scaling [37] and give
rise to a bound state BG. The addition of a small noise to the initial
condition(IC) facilitates the “phase mixing” of the finite-gap semiclasical
$\mathop{\rm dn}\nolimits$ potential, giving rise, in the long-time limit, to
a spatially uniform and statistically stationary integrable turbulence,
associated with a BG. The fundamental property of long-time “thermalization”,
or relaxation to a statistically stationary state, in an integrable system was
established numerically in several scenarios of the evolution of random waves
in the fNLS equation [40, 41, 42]. Here we extend this result to the
qualitatively new framework of the BG fission, which encompasses a broad range
of scenarios of transition to stationary integrable turbulence depending on
the value of $m$ in the semiclassical elliptic potential: from the nonlinear
development of the spontaneous modulational instability for $m\to 0$ to the
rarefied soliton gas fission for $m\to 1$. When the noise is sufficiently
small, the ZS spectrum remains essentially unchanged from that of the
semiclassical elliptic potential. The isospectrality of the BG fission enables
one to take advantage of the results of [21] and evaluate statistical measures
of the BG, such as the mean intensity and kurtosis, in terms of the spectral
density of states of the initial elliptic potential. In [21], it was predicted
that the kurtosis doubles in the long-time fNLS fission of the so-called
partially coherent waves into a SG. Here we show analytically that this result
generalizes to the BG fission from semiclassical elliptic potentials, and we
confirm our results by comparison with direct numerical simulations of the
time evolution of the noise-augmented potential.

An important ingredient of our construction is the interpretation of a BG as a
“composite SG” comprising of two spectrally distinct components: a regular SG
and a soliton condensate, defined as a critically dense SG [37] and describing
the modulationally unstable background in a BG. This composite SG provides a
natural extension of the “pure” fNLS soliton condensate framework used in [33,
43] to model the development of spontaneous (noise-induced) modulational
instability of fNLS plane wave and genus one elliptic solutions.

### NLS and semiclassical elliptic potentials

The cubic fNLS equation is, in normalized and dimensionless form and in the
semiclassical scaling,

$\displaystyle i\epsilon q_{t}+\epsilon^{2}q_{xx}+2\,|q|^{2}q=0\,,$ (1)

where subscripts $x$ and $t$ denote partial differentiation, $q(x,t)$
describes a complex-valued envelope of oscillations, and the physical meaning
of the variables $x$ and $t$ varies depending on the physical context. (E.g.,
in nonlinear fiber optics, $x$ is a retarded time and $t$ is the propagation
distance through the medium.) Equation (1) is the compatibility condition
${\mathbf{v}}_{xt}={\mathbf{v}}_{tx}$ of its Lax pair, namely the
overdetermined linear system

$\epsilon{\mathbf{v}}_{x}=(-iz\sigma_{3}+Q)\,{\mathbf{v}}\,,\quad\epsilon{\mathbf{v}}_{t}=P\,{\mathbf{v}}\,,$
(2)

for ${\mathbf{v}}(x,t,z)=(v_{1},v_{2})^{T}$, with
$P=-2iz^{2}\sigma_{3}+i(|q|^{2}+Q_{x})\,\sigma_{3}-2zQ$, where
$\sigma_{3}=\mathop{\rm diag}\nolimits(1,-1)$ is the third Pauli matrix, and

$Q(x,t)=\begin{pmatrix}0&q\\\ -q^{*}&0\end{pmatrix}\vspace*{-1ex}$ (3)

the asterisk denoting complex conjugate. The first half of Eq. (2) is the ZS
scattering problem, and $q$ and $z$ are referred to respectively as the
potential and the scattering parameter, or eigenvalue, since the first half of
Eq. (2) can be rewritten as the eigenvalue problem
${\mathcal{L}}{\mathbf{v}}=z{\mathbf{v}}$, where
${\mathcal{L}}:=i\sigma_{3}(\epsilon\partial_{x}-Q)$ is a one-dimensional
Dirac operator.

In [38] we characterized the Lax spectrum $\Sigma({\mathcal{L}})$ [44] of the
above ZS problem with the potential

$q(x,0)=\mathop{\rm dn}\nolimits(x;m)\,,\quad x\in\mathbb{R}\,,$ (4)

where $\mathop{\rm dn}\nolimits(x;m)$ is a Jacobi elliptic function [46, 45]
and $m\in(0,1)$ is the elliptic parameter. Note the spatial period is
$2L=2K_{m}$, where $K_{m}=K(m)$ is the complete elliptic integral of first
kind. Two limiting cases are $m=0,1$. The case $m=0$ reduces to the constant
potential $q(x,0)\equiv 1$, which is solved trivially, and for which
$\Sigma({\mathcal{L}})=\mathbb{R}\cup[-i,i]$. Conversely, as $m\to 1$,
$q(x,0)\to\mathop{\rm sech}\nolimits x$ and $K_{m}\to\infty$. This problem was
studied in [47], where it was shown that $\Sigma({\mathcal{L}})$ comprises the
real $z$-axis plus a set of purely imaginary discrete eigenvalues uniformly
distributed in the interval $(-i,i)$. Moreover, it was shown that $\mathop{\rm
sech}\nolimits x$ is a reflectionless potential (giving rise to pure soliton
solutions) if and only if $\epsilon=1/N$, with $N\in\mathbb{N}$. This Letter
connects these two limiting cases by characterizing the class of periodic
potentials (4) for all $0<m<1$. Figure 1 shows for $\epsilon=1/20$ and $m=0.7$
the time evolution of the dn potential according to fNLS without noise (left),
and with Gaussian noise (right).

Figure 1: Density plots of the time evolution of the potential in Eq. (4)
according to fNLS equation (1) with $\epsilon=1/20$ and: (a) $m=0.7$ without
noise; (b) $m=0.7$ with noise (see text for details). In both cases, the
spatial window coincides with the period of the IC.

### Band-gap structure

The structure of the Lax spectrum can be studied via Bloch-Floquet theory
using the monodromy matrix $M(z)=\Phi^{-1}(x,z)\Phi(x+2L,z)$, where
$\Phi(x,z)$ is a fundamental matrix solution of the scattering problem, and
$2L$ is the minimal period of the potential [48]. In turn, the Floquet
discriminant is defined by $\Delta(z)=\mathop{\rm tr}\nolimits M(z)/2$.
Bounded solutions of the scattering problem then exist for $z\in\mathbb{C}$
such that $\mathop{\rm Im}\nolimits\Delta(z)=0$ and $-1\leq\mathop{\rm
Re}\nolimits\Delta(z)\leq 1$. Importantly, $\Delta(z)$ is an entire function
of $z$ that is also Schwarz-symmetric: $\Delta(z^{*})=\Delta^{*}(z)$. As a
result, it is sufficient to consider the upper half plane, $\mathop{\rm
Im}\nolimits z\geq 0$. Since ${\mathcal{L}}$ is non-self-adjoint,
$\Sigma({\mathcal{L}})$ is not confined to the real $z$-axis. Nonetheless, the
zero-level curves of $\mathop{\rm Im}\nolimits\Delta(z)$ define a countable
set of analytic arcs $\Gamma_{n}$. Along each arc, the requirement
$-1\leq\mathop{\rm Re}\nolimits\Delta(z)\leq 1$ then defines a spectral band.
With these definitions, one can talk about bands as in a self-adjoint problem,
the difference being that in the non-self-adjoint case the bands are not
restricted to lie along the real $z$-axis, but lie instead along the
$\Gamma_{n}$.

The endpoints of the bands are periodic and antiperiodic eigenvalues of
${\mathcal{L}}$, namely, the values $z$ for which $\Delta(z)=\pm 1$, which
give rise to periodic and antiperiodic eigenfunctions, respectively. Moreover,
the real $z$-axis is an infinitely long band, and any band that intersects the
real $z$-axis transversally is called a “spine”. Any potential whose spectrum
is comprised of a finite number of bands is called a ”finite-band” potential
[8, 9, 49]. These special potentials correspond to the finite-genus solutions
of all equations in the Ablowitz-Kaup-Newell-Segur (AKNS) hierarchy.

### Elliptic finite-band potentials

In [38] it was proved that for the dn potential (4) and $\epsilon=1/N$,
$N\in\mathbb{N}$, the Lax spectrum of the ZS scattering problem comprises $2N$
Schwarz-symmetric bands along the interval $(-i,i)$, and produces a genus
$2N-1$ solution of fNLS [50]. Figure 2 (left) illustrates the dependence of
the spectrum on the elliptic parameter $m$ (horizontal axis), in particular,
the periodic (red curves) and antiperiodic (blue curves) eigenvalues of
${\mathcal{L}}$ along the imaginary $z$-axis (vertical axis) for
$\epsilon=1/7$, computed using finite truncations of Eqs. (A12). Notice all
gaps are closed when $m=0$ and they open as $m>0$ and remain open for all
$m\in(0,1)$. In the singular limit $m\to 1$, the band widths tend to zero, and
the periodic and antiperiodic eigenvalues “collide” to form the point spectrum
of ${\mathcal{L}}$ on the line. The corresponding spectrum in the complex $z$
plane for $m=0.9$ is shown in Fig. 2 (right).

Figure 2: Left: Periodic (red) and antiperiodic (blue) eigenvalues (vertical
axis) as a function of the elliptic parameter $m$ (horizontal axis) for
$\epsilon=1/7$. For every fixed $m\in(0,1)$, as $\epsilon\to 0$, the number of
bands grows like $1/\epsilon$. In the gray range, corresponding to
$\Gamma_{g}$, the band widths decrease exponentially in $\epsilon$.
Conversely, in the white range, corresponding to $\Gamma_{o}$, the band widths
decrease algebraically. Right: The spectrum of the dn potential (4) in the
complex $z$ plane with $m=0.9$ for $\epsilon=1/7$.

### Semiclassical limit and breather gas

Even though ${\mathcal{L}}$ is not self-adjoint, its spectrum in the limit
$\epsilon\to 0$ can be effectively characterized using WKB methods [51, 52],
which allows one to obtain precise asymptotic estimates for the location of
band edges, band widths, and gap widths [50]. In turn, these calculations show
that the semiclassical limit of the spectrum of the elliptic potential (4) is
compatible with the thermodynamic spectral scaling of a SG or a BG [37, 53].
All bands lie along a finite Schwarz-symmetric arc
$\Gamma=(-iq_{\mathrm{max}},iq_{\mathrm{max}})$ in the complex $z$-plane. As
$\epsilon\to 0$, however, the behavior of bands and gaps in the range
$\Gamma_{o}=(0,iq_{\mathrm{min}})$ differs from those in
$\Gamma_{g}=(iq_{\mathrm{min}},iq_{\mathrm{max}})$. Inside $\Gamma_{g}$, band
widths decay exponentially in $\epsilon$, the ratio of band width to gap width
tends to zero, and the bands accumulate according to some limiting density
$\varphi(z)$ [50]. Conversely, inside $\Gamma_{o}$, it is the relative band
gaps that tend to zero, giving rise to a continuous band in the limit. Thus,
exponentially small bands outside the finite band give rise to an “effective
soliton ensemble” on an “effective nonzero background”, i.e., a generalized
BG. As $m\to 1$, the range $\Gamma_{o}$ tends to zero, and one gets a pure SG.
Conversely, as $m\to 0$, it is the range $\Gamma_{g}$ that tends to zero, and
one obtains a pure condensate. We therefore have that, in the semiclassical
limit, the potential (4) gives rise to a composite SG. Next we characterize
this gas in detail.

We can derive the two key quantities in the spectral theory of soliton and
breather gases, namely, the density of bands, and the normalized logarithmic
band width. In particular, the density of bands $\varphi(z)$ is [50]

$\varphi(z)=\frac{2|z|}{\pi}\int_{0}^{x_{o}(z)}[\mathop{\rm
dn}\nolimits^{2}(x;m)-|z|^{2}]^{-1/2}\,\mathrm{d}x\,,$ (5)

where $x_{o}(z)=K_{m}$ for $z\in(0,iq_{\mathrm{min}})$ and $x_{o}(z)>0$ is a
simple turning point [at which $\mathop{\rm dn}\nolimits(x;m)=|z|$] for
$z\in(iq_{\mathrm{min}},i)$ [54]. Following [37, 23, 53], we obtain the
density of states (DOS) $f(z)$ as [55, 50] $f(z)=\varphi(z)/(2K_{m})$. As
$m\to 0$, $\Gamma_{g}$ collapses, and $f(\eta)$ tends to the Weyl distribution
of a soliton condensate [50]. Conversely, as $m\to 1$, $L\to\infty$,
$\Gamma_{o}$ collapses, and $f(\eta)$ tends to zero, as appropriate for a
rarified SG. The interpretation of a BG as a composite soliton gas is also
supported by the results of [59], where basic fNLS breathers were generated
numerically with high accuracy by appropriately configuring $N$-soliton fNLS
solutions with large $N$.

### Kurtosis

The DOS allows one to obtain the following expression for the kurtosis
$\kappa$. In [21], the following general expression for the kurtosis $\kappa$
of spatially uniform bound state integrable turbulence was obtained:

$\kappa=\frac{\langle|q|^{4}\rangle}{\langle|q|^{2}\rangle^{2}}=\frac{2\overline{\eta^{3}}}{3\overline{\eta}^{2}}\,,$
(6)

where the angle brackets denote ensemble average, $\eta=\mathop{\rm
Im}\nolimits z$, and the overbar denotes average over the DOS:

$\overline{\eta^{k}}=\int_{\gamma}\eta^{k}f(\eta)\,\mathrm{d}\eta=(-1)^{\frac{k+1}{2}}\frac{c_{k}}{L}\int_{0}^{L}q^{k+1}(x,0)\,\mathrm{d}x\,$
(7)

are the moments of the DOS $f(\eta)$, where $\gamma$ is the spectral support
of the SG, and where the last equality follows from [53], with
$c_{1}=\frac{1}{4}$ and $c_{3}=\frac{3}{16}$. Thus,

$\displaystyle\kappa=2\kappa_{0},$ (8a)
$\displaystyle\kappa_{0}=L\frac{\int_{0}^{L}|q(x,0)|^{4}\,\mathrm{d}x}{\big{(}\int_{0}^{L}|q(x,0)|^{2}\,\mathrm{d}x\big{)}^{2}}\,.$
(8b)

This result is consistent with the kurtosis doubling obtained in [21] for the
SG fission from partially coherent waves, with a crucial difference: here
$\kappa_{0}$ is computed from a purely deterministic IC. Moreover, the above
results imply that [50] any periodic SG or BG of the fNLS equation generated
by a deterministic real and even single-lobe potential $q(x,0)$ has kurtosis
$\kappa\geq 2$, with $\kappa=2$ if and only if $q(x,0)$ is constant. Recall
that $\kappa>2$ implies heavy-tailed non-Gaussian statistics, which is an
indication of the possible presence of rogue waves [56, 57, 58]. In
particular, for the potential (4), straightforward calculations yield

$\kappa_{0}=K_{m}[2(2-m)E_{m}-(1-m)K_{m}]/(3E_{m}^{2})\,,$ (9)

with $E_{m}$ the complete elliptic integral of the second kind.

### Temporal dynamics and numerical validation

Next we discuss the numerically computed time evolution of the fNLS solutions
produced by the potential (4). Figure 1 demonstrates how the addition of a
small amount of zero-mean white Gaussian noise (with standard deviation
$\sigma=10^{-2}$) to the IC effectively randomizes the solution. Figure
3(left) shows the temporal evolution of the normalized fourth moment
$|q|^{4}/(|q|^{2})^{2}$ of the resulting fNLS solutions as well as their
ensemble average (i.e., the kurtosis), clearly demonstrating the establishment
of a statistically stationary integrable turbulence state. Finally Fig.
3(right) shows how the dependence of the kurtosis on the elliptic parameter
$m$ matches the theoretical prediction from Eqs. (6–9), providing a strong
validation of the analytical predictions.

Figure 3:  Left: The averaged temporal evolution of the kurtosis from an
ensemble of noise realizations (red line). The gray region shows the envelope
of the normalized fourth moment of the solution from all individual noise
realizations. Right: The kurtosis of a fully randomized gas as a function of
the elliptic parameter $m$. Blue curve: theoretical prediction from Eqs. (8a)
and (9). Dashed black curve: $\kappa_{0}$ from Eq. (9). Red squares:
Numerically computed asymptotic value of kurtosis. Yellow stars: Ensemble
average of the numerically evaluated integral from Eq. (8b) for ICs consisting
of the dn potential plus noise.

### Concluding remarks

In summary, we have presented an analytically tractable model describing a
mechanism for the formation of integrable turbulence via BG fission of a
semiclassical elliptic potential augmented by weak noise. Our analysis
introduces a natural interpretation of a BG as a “composite SG”, consisting of
two distinct components: a regular SG plus a soliton condensate. The
analytical model comprises a one-parameter family of such BGs, which
interpolates between a pure soliton condensate (as $m\to 0$) and a rarified SG
(as $m\to 1$) [60, 61, 62, 63]. Intermediate values of $m$ give rise to a
mixed regime that interpolates between the above two extremes. We validated
the theoretical results with direct numerical simulations of the BG fission
and demonstrated the establishment, at large $t$ of a statistically stationary
integrable turbulence field characterised by the kurtosis value which is shown
analytically and confirmed numerically to be twice as large as the fourth
normalized moment of the initial elliptic potential, implying, in particular,
the presence of rogue waves for all nonzero values of $m$.

We expect that similar results will hold for arbitrary real single-lobe
periodic potentials. Our results also open up a number of interesting avenues
for further research. Since the ZS problem is common to all equations of the
AKNS hierarchy, the results of this Letter will generate SGs/BGs for all such
equations, which although spectrally equivalent, exhibit qualitatively
different dynamics. The study of their resulting SGs/BGs is threfore an
interesting open question. For example, despite some recent work [64, 65], the
spectral theory of SGs for the focusing modified Korteweg-deVries (mKdV)
equation is still open, and we expect that the mKdV soliton gas phenomenology
will be different than that of both the KdV and focusing NLS gases,
potentially involving solitons of both polarities in the same gas.

### Acknowledgement

We acknowledge fruitful discussions at the Dispersive Hydrodynamics Program
hosted by the Issac Newton Institute for Mathematical Sciences in 2022. We
also thank Mark Hoefer for useful suggestions. The work of G.B. and A.T. was
partially supported by the National Science Foundation under grants
DMS-2009487 and DMS-2009647, respectively, and that of X.D.L. by the National
Natural Science Foundation of China under grant number 12101590.

## Appendix

In this appendix we provide additional details on the results presented in the
main text.

### Jacobi dn potential and the Lax spectrum

In [38] we studied the focusing NLS equation

$iq_{t}+q_{xx}+2|q|^{2}q=0,$ (A1)

where the initial condition (IC) is a multiple Jacobi “dn” elliptic function
[45]:

$q(x,0)=A\mathop{\rm dn}\nolimits(x;m),\quad A\in\mathbb{R}.$ (A2)

As $m$ goes from 0 to 1 the fNLS dynamics corresponding to $N\mathop{\rm
dn}\nolimits(x;m)$ initial data interpolates between plane-wave background
($m=0$), to genus $2N-1$ finite-gap solution ($0<m<1$), to pure $N$-soliton
($m=1$).

To connect the above problem to the semiclassical setting, it is sufficient to
rescale $q$ and $t$ to make the initial data independent of $A$, by letting
$q(x,t)\mapsto Aq(x,At)$, which yields

$\displaystyle i\epsilon q_{t}+\epsilon^{2}q_{xx}+2|q|^{2}q=0,$ (A3)
$\displaystyle q(x,0)=\mathop{\rm dn}\nolimits(x;m),$ (A4)

with $\epsilon=1/A$. Hence, studying the semiclassical limit ($\epsilon\to 0$)
of Eq. (A3) with IC (A4) is equivalent to studying the large-$A$ limit of Eqs.
(A1)–(A2).

Figure A2(left) depicts the potential (A4) for various values of the elliptic
parameter $m\in(0,1)$. Recall that the real period of the dn potential is
$2K_{m}$, where $K_{m}=K(m)$ is the complete elliptic integral of the first
kind [45]. Figure A2(right) shows the contours $\\{z:\mathop{\rm
Im}\nolimits\Delta(z)=0\\}$ (thin black curves) together with the Lax spectrum
(thick blue curves) for a generic (non-integer) value of the semiclassical
parameter $\epsilon$. Figure A2 shows the noise-augmented initial conditions
(ICs) for two different values of $m$ and Fig. A3 compares the Lax spectrum
with and without noise, demonstrating that the noise does not affect it
appreciably. In particular, the number and locations of the solitonic
excitations remains unchanged.

In [38] we proved that the Lax spectrum of the focusing ZS scattering problem
associated with Eq. (A1) and dnoidal potential (A2) and $A\in\mathbb{Z}$ is
comprised of $2A$ Schwarz-symmetric bands along the imaginary axis of the
spectral variable. As a result, the corresponding solutions of focusing NLS
are special “finite-gap” solutions of genus $2A-1$. Next, we provide some
details on how the spectrum can be efficiently computed numerically.

Figure A1: Left: $\mathop{\rm dn}\nolimits(x;m)$ (vertical axis) as a function
of $x$ (horizontal axis) for various values of the elliptic parameter $m$.
From top to bottom: $m=0$ (constant), 0.6 (gray), 0.9 (black), 0.99 (gray) and
0.999 (red). Right: Diagram illustrating the curves $\mathop{\rm
Im}\nolimits\Delta(z)=0$ (black) and the Lax spectrum $\Sigma({\mathcal{L}})$
(blue) in the complex $z$ plane for the potential (A4) with generic (non-
integer) value of $\epsilon$.

Figure A2: Red curve: Initial condition (A4). Black curve: Same augmented by
Gaussian noise with mean $\mu=0$ and standard deviation $\sigma=10^{-2}$.
Left: $m=0.05$. Right: $m=0.9$.

### PT-symmetric periodic Schrödinger potentials with real spectrum

When the potential is real-valued as with Eq. (A4), $\Delta(z)$ also possesses
an additional symmetry: $\Delta(-z^{*})=\Delta(z)$, and thus the eigenvalues
come in quartets, i.e., $\\{z,z^{*},-z,-z^{*}\\}$. Moreover, the
transformation ${\mathbf{w}}=\Lambda{\mathbf{v}}$, where
${\mathbf{w}}(x,\lambda)=(w_{+},w_{-})^{T}$ and

$\Lambda=\frac{1}{\sqrt{2}}\begin{pmatrix}[r]1&i\\\ 1&-i\end{pmatrix}\,,$ (A5)

maps the ZS problem into the time-independent Schrödinger equation with a
complex periodic potential, namely,

$H_{\pm}w_{\pm}=\lambda w_{\pm}\,,\quad\lambda=z^{2}\,,$ (A6a) where
$H_{\pm}=-\epsilon^{2}\partial_{x}^{2}+V_{\pm}(x)\,,\quad V_{\pm}(x)=\mp
i\epsilon q_{x}-q^{2}\,.$ (A6b)

Recall that the Schrödinger equation with a periodic potential is referred to
as Hill’s equation [66]. If $q$ is real and even, then the potentials are PT-
symmetric: $V_{\pm}(-x)=V_{\pm}^{*}(x)$. In particular, for the elliptic
potential Eq. (A4), $H_{\pm}$ is

$H_{\pm}=-\epsilon^{2}\partial_{x}^{2}-\mathop{\rm dn}\nolimits^{2}(x)\pm
i\epsilon m\mathop{\rm sn}\nolimits(x)\mathop{\rm cn}\nolimits(x)\,.$ (A7)

Since $\Sigma(H_{+})=\Sigma(H_{-})$ it is enough to consider only $H_{-}$.
Below we discuss that, even though $H_{\pm}$ are not Hermitian, all of their
eigenvalues are real.

A further simplification is obtained via the change of variable
$y=2\mathop{\rm am}\nolimits(x;m)$ [where $\mathop{\rm am}\nolimits$ is the
Jacobi amplitude], which maps Eq. (A7) into a complex perturbation of Ince’s
equation [66]:

$4\epsilon^{2}[1-m\sin^{2}(y/2)]\,w_{yy}-\epsilon^{2}(m\sin y)\,w_{y}\\\
+[\lambda+(1-m\sin^{2}(y/2))+i\epsilon{\textstyle\frac{m}{2}}\,\sin
y]\,w=0\,.$ (A8)

Bloch-Floquet theory implies that any bounded solutions of Eq. (A8) can be
written as $w(y)=\mathop{\rm e}\nolimits^{i\nu y}p(y)$, where
$\nu\in\mathbb{R}$ is the Floquet exponent and $p(y+2\pi)=p(y)$, and can
therefore be expanded in Fourier series, with Fourier coefficients given by a
three-term recurrence relation. The eigenvalues of the ZS problem are then
related to the (real or) complex values of $\lambda$ for which the above ODE
admits bounded solutions. Integer and half-integer values of $\nu$ yield
respectively periodic and antiperiodic eigenfunctions.

Figure A3: Left: Numerically computed spectrum of the scattering problem for
the dn potential with $m=0.9$ via the Fourier-Hill method [70]. Right: Same
but with noise added to the potential with $\mu=0$ and $\sigma=10^{-2}$.

### Three-term recurrence relation

Any solution of Eq. (A8) that is bounded for all $y\in\mathbb{R}$ can be
represented by a Fourier series as

$w(y;\lambda)=\mathop{\rm e}\nolimits^{i\nu
y}\sum_{n\in\mathbb{Z}}c_{n}\mathop{\rm e}\nolimits^{iny}\,,$ (A9)

with $\nu\in\mathbb{R}$. The coefficients $\\{c_{n}\\}_{n\in\mathbb{Z}}$ are
given by the three-term recurrence relation

$\alpha_{n}c_{n-1}+(\beta_{n}-\lambda)c_{n}+\gamma_{n}c_{n+1}=0\,,$ (A10)

where

$\displaystyle\alpha_{n}=-m\,[{\textstyle\frac{1}{2}}-\epsilon(n+\nu-1)][{\textstyle\frac{1}{2}}+\epsilon(n+\nu-{\textstyle\frac{1}{2}})]\,,$
(A11a) $\displaystyle\beta_{n}=(1-{m}/{2})[\epsilon^{2}(2n+2\nu)^{2}-1]\,,$
(A11b)
$\displaystyle\gamma_{n}=-m\,[{\textstyle\frac{1}{2}}-\epsilon(n+\nu+1)][{\textstyle\frac{1}{2}}+\epsilon(n+\nu+{\textstyle\frac{1}{2}})]\,,$
(A11c)

$n\in\mathbb{Z}$.

### Elliptic finite-band potentials

When $\epsilon=1/N$, some of the coefficents in Eqs. (A11) vanish, and as a
result it is possible to decompose the doubly-infinite recurrence relation
into two semi-infinite ones. Specifically, let

$\displaystyle B_{o}^{\pm}$
$\displaystyle=\begin{pmatrix}\beta_{j}&\gamma_{j}\\\
\alpha_{j+1}&\beta_{j+1}&\gamma_{j+1}\\\ &\ddots&\ddots&\ddots\end{pmatrix},$
(A12a) $\displaystyle B_{\infty}^{\pm}$
$\displaystyle=\begin{pmatrix}\beta_{-j}&\alpha_{-j}\\\
\gamma_{-j-1}&\beta_{-j-1}&\alpha_{-j-1}\\\
&\ddots&\ddots&\ddots\end{pmatrix},$ (A12b)

with $\nu=N/2$ and $j=0$ for the minus sign and $\nu=(1-N)/2$ and $j=1$ for
the plus sign. The periodic and antiperiodic spectrum of $H_{-}$ is the union
of the spectra of $B_{o}^{\pm}$ and $B_{\infty}^{\pm}$. Specifically: when $N$
is even, $B_{o}^{-}$ and $B_{\infty}^{-}$ yield the periodic eigenvalues and
$B_{o}^{+}$ and $B_{\infty}^{+}$ the antiperiodic ones, and viceversa when $N$
is odd. In [38] we proved that the eigenvalues of all four of these half-
infinite matrices are real. This is the key to prove that, for any
$m\in(0,1)$, the potential $q$ in Eq. (A4) for the focusing ZS scattering
problem is finite-band if and only if $\epsilon=1/N$ with $N\in\mathbb{N}$.
(The result is easily extended to $N\in\mathbb{Z}$ by phase invariance.)
Moreover, if $\epsilon=1/N$, $q$ is a $2N$-band potential, and

$\Sigma({\mathcal{L}})=\mathbb{R}\cup\Big{(}\\!\\!\bigcup_{n=1}^{N}[-i\eta_{2n},-i\eta_{2n-1}]\cup[i\eta_{2n-1},i\eta_{2n}]\Big{)}\,,$
(A13)

where $0<\eta_{1}<\eta_{2}<\cdots<\eta_{2N}<1$. This implies that the spectral
curve and the flow induced by each member of the Ablowitz-Kaup-Newell-Segur
(AKNS) hierarchy has finite genus $2N-1$. Importantly, the finite truncations
of the matrices (A12) provide an efficient way to numerically compute the
spectrum.

### Semiclassical WKB analysis of the ZS problem

Next we consider the semiclassical limit of the spectrum, namely the limiting
behavior of bands and gaps as $\epsilon\to 0$. It is useful to briefly recall
the asymptotic analysis of the focusing ZS scattering problem via WKB methods
from [51]. We emphasize that the results of this and the next few sections
apply to a broad class of potentials, not just dn.

Suppose that $q(x)$ is the $2L$-periodic extension of a real, even, non-
negative single-lobe potential. Thus, $q$ has one maximum and one minimum in
$(-L,L]$, which without loss of generality can be taken to be respectively at
$x=0$ and $x=L$. To avoid trivial cases, assume that $q$ is not constant. Let

$q_{\rm max}=q(0)\,,\quad q_{\rm min}=q(L)\,,$ (A14a) and
$s(x,z)=\int_{-x_{o}(z)}^{x}\sqrt{|q^{2}(u)+z^{2}|}\,\mathrm{d}u\,,$ (A14b)

where $z\in\mathbb{R}\cup i\mathbb{R}$, and $x_{o}(z)$ is a simple (real)
turning point. Without loss of generality we can limit ourselves to
considering $\mathop{\rm Im}\nolimits z\geq 0$ thanks to the Schwarz symmetry
of the spectrum. Then, as $\epsilon\to 0$ (see [51, 52]):

  * (i)

For $z\in\Gamma_{o}\cup\mathbb{R}$, where $\Gamma_{o}=(0,iq_{\rm min})$, one
has $q^{2}(x)+z^{2}>0$ for $x\in(-L,L)$. There are no turning points, and

$\Delta(z;\epsilon)\sim\cos(s_{\rm i}(z)/\epsilon)\,,$ (A15)

where $s_{\rm i}(z)=\int_{-L}^{L}\sqrt{q^{2}(x)+z^{2}}\,\mathrm{d}x$. This
range corresponds to the region in white in Fig. 2.

  * (ii)

For $z\in\Gamma_{g}$, where $\Gamma_{g}=(iq_{\rm min},iq_{\rm max})$, there
are two real symmetric turning points $\pm x_{o}$, i.e., values at which
$q^{2}(\pm x_{o})+z^{2}=0$. In this region (see Fig. A4)

$\Delta(z;\epsilon)\sim\cos(s_{1}(z)/\epsilon)\cosh(s_{2,\epsilon}(z)/\epsilon)\,,$
(A16)

where $s_{1}(z)=s(x_{o},z)$ and $s_{2,\epsilon}(z)=\epsilon\ln 2-2s(-L,z)$.
This range corresponds to the region in gray in Fig. 2.

  * (iii)

For $z\in(iq_{\rm max},i\infty)$ one has $q^{2}(x)+z^{2}<0$ for $x\in(-L,L)$.
There are no turning points and

$\Delta(z;\epsilon)\sim\cosh(s_{\rm ii}(z)/\epsilon)\,,$ (A17)

where $s_{\rm ii}(z)=\int_{-L}^{L}\sqrt{-q^{2}(x)-z^{2}}\,\mathrm{d}x$.

Figure A4 shows a comparison between the predicted behavior of
$\Delta(z,\epsilon)$ and the numerically computed value.

Figure A4: Left: the numerically computed Floquet discriminant (in blue) for
the potential (A4) with $m=0.96$ and $\epsilon=1/20$, using fourth-order
Runge-Kutta versus the WKB approximation (in red). The black solid horizontal
lines show $\Delta=\pm 1$. The black dashed vertical lines show $\mathop{\rm
Im}\nolimits(z)=q_{\rm min}$ and $\mathop{\rm Im}\nolimits(z)=q_{\rm max}$.
Note that, to scale the plot, the function $F(\tau)=\tau$ for $|\tau|\leq 1$
and $F(\tau)={\rm sgn}(\tau)(1+\log_{10}|\tau|)$ for $|\tau|>1$ is plotted.
Right: The number of spectral bands for $z\in\Gamma_{g}$ as predicted by WKB
(red curve) versus the numerically computed values (blue triangles) as a
function of $\epsilon$.

### Semiclassical limit of bands and gaps in $\bm{\Gamma}_{g}$

When $z\in(iq_{\rm min},iq_{\rm max})$, Eq. (A16) implies that the bands are
approximately centered at the roots $\\{z_{n}\\}$ of the equation
$s_{1}(z_{n})=(n-{\textstyle\frac{1}{2}})\pi\epsilon$. Thus, the asymptotic
number of bands for $z\in(iq_{\rm min},iq_{\rm max})$ is

$N^{\epsilon}_{\mathrm{bands}}\sim\Big{\lfloor}\frac{J_{g}}{\pi\epsilon}\Big{\rfloor},\qquad\epsilon\to
0\,,$ (A18a) where $J_{g}=s_{1}(iq_{\rm
min})=\int_{-L}^{L}\sqrt{q^{2}(x)-q^{2}_{\rm min}}\,\mathrm{d}x\,.$ (A18b)

and $\lfloor\cdot\rfloor$ is the floor function.

Next, consider the limit $\epsilon\to 0$ and $n\to\infty$ satisfying
$0<n\epsilon<J_{g}/\pi$ and $n$ is the band index. Denote the $n$-th band
width as $w_{n}^{\epsilon}=|z_{n}^{+}-z_{n}^{-}|$, where
$\Delta(z_{n}^{\pm})=\pm 1$, respectively. Similarly, denote the $n$-th gap
width as $g_{n}^{\epsilon}=|z_{n+1}^{\pm}-z_{n}^{\pm}|$. Using
$\Delta(z_{n}^{\pm})=\pm 1$ together with Eq. (A16), one gets

$w_{n}^{\epsilon}\sim\frac{4\epsilon}{|s_{1}^{\prime}(z_{n})|}\mathop{\rm
e}\nolimits^{-s_{2,\epsilon}(z_{n})/\epsilon},\qquad\epsilon\to 0\,.$ (A19)

Moreover, since $g_{n}^{\epsilon}\sim|z_{n+1}-z_{n}|$ and
$s_{1}(z_{n+1})-s_{1}(z_{n})=\pi\epsilon$ it follows that

$g_{n}^{\epsilon}\sim\frac{\pi\epsilon}{|s_{1}^{\prime}(\xi_{n})|},\qquad\epsilon\to
0\,,$ (A20)

where $\mathop{\rm Im}\nolimits(z_{n+1})<\mathop{\rm
Im}\nolimits(\xi_{n})<\mathop{\rm Im}\nolimits(z_{n})$. Finally, the band-to-
gap ratio $w_{n}^{\epsilon}/g_{n}^{\epsilon}\to 0$ as $\epsilon\to 0$
exponentially fast in this region.

### Effective solitons in the fNLS equation

Recall that, for localized potentials (i.e., $q(x,t)\to 0$ rapidly as
$|x|\to\infty$), solitons are parameterized by the discrete eigenvalues of the
ZS scattering problem. In the periodic problem, nonlinear excitations are
considered to be “effective solitons” if the relative band width is less than
some small threshold parameter [51, 68, 67], i.e., if $W_{n}<\kappa$, where
$0<\kappa\ll 1$ and

$W_{n}^{\epsilon}=\frac{w_{n}^{\epsilon}}{w_{n}^{\epsilon}+g_{n}^{\epsilon}}\,.$
(A21)

Thus, rewriting the relative band width as

$W_{n}^{\epsilon}=\frac{w_{n}^{\epsilon}}{g_{n}^{\epsilon}}\Big{(}\frac{1}{1+w_{n}^{\epsilon}/g_{n}^{\epsilon}}\Big{)}\,,$
(A22)

one can show that

$W_{n}^{\epsilon}\sim\frac{4}{\pi}\mathop{\rm
e}\nolimits^{-s_{2,\epsilon}(z_{n})/\epsilon},\qquad\epsilon\to 0\,.$ (A23)

Note that for $z\in(iq_{\rm min},iq_{\rm max})$ the relative band width
$W_{n}$ is monotonic decreasing. Thus, the effective solitons are confined to
the interval

$z\in(z_{s},iq_{\rm max})\subset(iq_{\rm min},iq_{\rm max})\,,$ (A24)

where $z_{s}$ is the unique solution to $W_{n}=\kappa$. For $z\in(iq_{\rm
min},iq_{\rm max})$ let $s_{2}(z)=-s(-L,z)\geq 0$ (see Eq. (A14b)), then to
leading order $z_{s}$ is given implicitly by

$s_{2}(z_{s})=\frac{\epsilon}{2}\ln\Big{(}\frac{2}{\pi\kappa}\Big{)}\,,$ (A25)

and the number of effective solitons is given by

$N_{\mathrm{solitons}}=\Big{\lfloor}\frac{s_{1}(z_{s})}{\pi\epsilon}+\frac{1}{2}\Big{\rfloor}\,.$
(A26)

Expanding $s_{2}(z)$ about $z=iq_{\rm min}$ and evaluating at $z_{s}$, for any
$\kappa>0$ it follows that $z_{s}\to iq_{\rm min}$ as $\epsilon\to 0$. Thus,
in the semiclassical limit the entire interval $(iq_{\rm min},iq_{\rm max})$
is comprised of an infinite ensemble of effective solitons.

### Semiclassical limit of bands and gaps in $\Gamma_{o}$

Recall that, when $q=\mathop{\rm dn}\nolimits(x;m)$ with $m\in(0,1)$ and
$\epsilon=1/N$ with $N\in\mathbb{N}$ there are precisely $N$ bands for
$z\in(0,iq_{\max})$. Hence, not all bands become effective solitons in the
limit $\epsilon\to 0$ (see Fig. 2). For $z\in\Gamma_{o}$, following similar
arguments as for $z\in\Gamma_{g}$, one obtains

$w_{n}^{\epsilon}=O(\epsilon),\qquad\epsilon\to 0\,.$ (A27)

Thus, at $z_{*}=iq_{\min}=i\mathop{\rm dn}\nolimits(K_{m};m)$, which is the
boundary between the intervals $\Gamma_{g}=(iq_{\min},iq_{\max})$ and
$\Gamma_{o}=(0,iq_{\min})$, there is a transition from exponentially decaying
band widths (giving rise to a soliton gas), to algebraically decaying band
widths (giving rise to a soliton condensate). Thus, we have the physically
realistic scenario of a generalized breather gas — namely, a soliton gas on a
condensate background.

The number of bands with support in $\Gamma_{o}$ can be computed by noting
that Eq. (A15) holds in this interval with exponentially small corrections at
$\Delta(z)=\pm 1$. Thus, bands in this region are approximately centered at
the roots $\\{z_{n}\\}$ of the equation $s_{\rm
i}(z_{n})=(n-{\textstyle\frac{1}{2}})\pi\epsilon$. Thus, the asymptotic number
of bands for $z\in(0,iq_{\rm min})$ is

$N_{\rm
bands}^{\epsilon}\sim\Big{\lfloor}\frac{J_{o}}{\pi\epsilon}\Big{\rfloor},\qquad\epsilon\to
0\,,$ (A28a) where $J_{o}=s_{\rm i}(0)-s_{1}(iq_{\rm min}),$ (A28b) and
$s_{\rm i}(0)=\int_{-L}^{L}q(x)\,\mathrm{d}x\,.$ (A28c)

### Spectral characteristics of the composite soliton gas

Recall that the density of bands $\varphi(z)$ and scaled logarithmic band
width $\nu(z)$ are two key quantities in the spectral theory of soliton gases.
Following [53] we have that, for $z\in\Gamma_{g}$, and the dn potential (A4)

$\displaystyle\varphi(z)$ $\displaystyle=\frac{2|z|}{J}\int_{0}^{\mathop{\rm
dn}\nolimits^{-1}(|z|)}\frac{\mathrm{d}x}{\sqrt{\mathop{\rm
dn}\nolimits^{2}(x;m)+z^{2}}}\,,$ (A29a) $\displaystyle\nu(z)$
$\displaystyle=\frac{2\pi}{J}\int_{\mathop{\rm
dn}\nolimits^{-1}(|z|)}^{K_{m}}\sqrt{|\mathop{\rm
dn}\nolimits^{2}(x;m)+z^{2}|}\,\mathrm{d}x\,,$ (A29b)

with $J$ is a suitable normalization constant, determined below. Thus we can
easily express the spectral scaling function as

$\sigma(z)=\frac{2\nu(z)}{\varphi(z)}.$ (A30)

Similarly, for $z\in\Gamma_{o}$, we have

$\displaystyle\varphi(z)$
$\displaystyle=\frac{2|z|}{J}\int_{0}^{K_{m}}\frac{\mathrm{d}x}{\sqrt{\mathop{\rm
dn}\nolimits^{2}(x;m)+z^{2}}}\,,$ (A31a) $\displaystyle\sigma(z)$
$\displaystyle=0\,.$ (A31b)

(In $\Gamma_{o}$, the spectral scaling function $\sigma(z)$ is zero as
expected for a condensate.) Finally, ensuring that the integral of the density
of bands over the whole support $(0,i)$ is one, one can then show that [53]

$J=J_{g}+J_{o}=\int_{-K_{m}}^{K_{m}}\mathop{\rm
dn}\nolimits(x;m)\,\mathrm{d}x=\pi\,.$ (A32)

### Limits $m\to 0$ and $m\to 1$

Next we discuss the density of bands formulae in the limits $m\to 0^{+}$, and
$m\to 1^{-}$.

Recall that $\mathop{\rm dn}\nolimits(x;0)\equiv 1$ and $K(0)=\pi/2$. Thus,
$q_{\rm min}\to 1$ as $m\to 0^{+}$ and Eq. (A31a) has support $z\in(0,i)$.
Moreover, $J_{o}\to\pi$ as $m\to 0^{+}$. Thus,

$\varphi(z)\to\frac{|z|}{\sqrt{1+z^{2}}},\qquad m\to 0^{+}\,,$ (A33)

which corresponds to Weyl’s distribution as expected (see [33]).

Next, recall that $\mathop{\rm dn}\nolimits(x;1)\equiv\mathop{\rm
sech}\nolimits(x)$ and $K(m)\to\infty$ as $m\to 1^{-}$. Thus, $q_{\rm min}\to
0$ as $m\to 1^{-}$ and Eq. (A29a) has support $z\in(0,i)$. Moreover,
$J_{g}\to\pi$ as $m\to 1^{-}$, and $\mathop{\rm dn}\nolimits^{-1}(|z|)\to{\rm
arcsech}(|z|)$ as $m\to 1^{-}$. Thus, in this limit we get a uniform
distribution:

$\varphi(z)\to 1,\qquad m\to 1^{-}\,.$ (A34)

### Nonlinear dispersion relations

Together with the set $\Gamma_{+}=\Gamma_{g}\cup\Gamma_{o}$, the spectral
scaling function $\sigma(z)$ determines the integral equation (the first
nonlinear dispersion relation, NDR) DOS $f(z)$ for soliton gases:

$\displaystyle\int_{\Gamma_{+}}\log\left|\frac{\zeta-z^{*}}{\zeta-z}\right|f(\zeta)\,\mathrm{d}|\zeta|+\sigma(z)f(z)=\mathop{\rm
Im}\nolimits z.$ (A35)

(Note that in the portion of the spectrum corresponding to a soliton
condensate, i.e., $\Gamma_{o}$, one has $\sigma(z)=0$.) The second NDR
involves the density of fluxes, but this quantity is zero in our case since
all the nonlinear excitations have zero velocity. It was proven in [53] that
$f(z)=r\varphi(z)$ solves the first NDR for soliton gases, where

$r=\frac{J}{2\pi L}=\frac{1}{2K_{m}}\,.$ (A36)

As a result, expressions (A29a) and (A31a) also yield the density of states
(DOS) $f(z)$ for the dn potential (A4).

### Kurtosis and Jensen’s inequality

In the main text,we have shown that the kurtosis $\kappa$ of a fully developed
SG generated by a deterministic real and even single-lobe initial condition
$q(x,0)$ is given by $\kappa=2\kappa_{0}$, where

$\kappa_{0}=\frac{L\int_{0}^{L}|q(x,0)|^{4}\,\mathrm{d}x}{\big{(}\int_{0}^{L}|q(x,0)|^{2}\,\mathrm{d}x\big{)}^{2}}\,.$
(A37)

Next, the rescaling $y=x/L$ and $\tilde{q}(y,0)=q(Ly,0)$ implies

$\kappa_{0}=\frac{\int_{0}^{1}|\tilde{q}(y,0)|^{4}\,\mathrm{d}y}{\big{(}\int_{0}^{1}|\tilde{q}(y,0)|^{2}\,\mathrm{d}y\big{)}^{2}}\,.$
(A38)

Then Jensen’s inequality [45] implies that $\kappa_{0}\geq 1$, with the
equality holding only if $q(x,0)$ is constant.

$m$ | $\kappa_{o}$ | $\kappa_{\infty}$ |  $T_{\infty}$ |  $~{}M$ |  $N_{\mathrm{Fourier}}$
---|---|---|---|---|---
$0.05$ | 1.007 | 1.92 | $4$ | 500 | 1024
$0.3$ | 1.016 | 1.96 | $4.5$ | 400 | 1024
$0.5$ | 1.059 | 2.07 | $4$ | 400 | 2048
$0.7$ | 1.17 | 2.37 | $2.5$ | 200 | 4096
$0.8$ | 1.29 | 2.66 | $2.5$ | 200 | 4096
$0.9$ | 1.53 | 3.022 | $3.0$ | 400 | 4096
Table 1: Numerical simulation parameters and results, where
$N_{\mathrm{Fourier}}$ is the number of Fourier modes used, $M$ is the number
of simulations (each of which corresponds to an independent noise
realization), and $T_{\infty}$ is the numerically observed approximate
thermalization time.

Figure A5: Time evolution of the kurtosis with $\epsilon=1/20$ and
$\sigma=10^{-2}$. Gray regions: Envelope of 100 ensemble realizations. Red
curves: Ensemble average of $M$ realizations. Top left: $m=0.05$ ($N=500$).
Top right: $m=0.3$ ($N=400$). Middle left: $m=0.5$ ($N=400$). Middle right:
$m=0.7$ ($N=200$). Bottom left: $m=0.8$ ($N=200$). Bottom right: $m=0.9$
($N=400$).

Figure A6: Numerical solution of fNLS with dn IC, with and without noise, and
$\epsilon=1/20$. Column 1: $\sigma=0$ (without noise); Column 2:
$\sigma=10^{-2}$ (with noise); Row 1: Density plot, $m=0.05$; Row 2: Solution
at $t=0.05$ and $t=4$ for $m=0.05$.

### Numerical methods

We performed several numerical experiments related to the kurtosis of randomly
generated solutions to the fNLS equation. In particular, we studied the time
evolution of fNLS in the semiclassical limit with elliptic IC plus a complex
perturbation of independent and identically distributed normal random
variables with mean $\mu=0$ and standard deviation $\sigma=10^{-2}$.

Due to modulational instability of fNLS the numerical computation of solutions
is a delicate topic, especially in the semiclassical limit. To this end, in
this work all numerical simulations of the fNLS equation were performed using
an eighth-order split-step method [69]. Moreover, to ensure numerical accuracy
the isospectral property of the scattering data was confirmed at several data
points as the simulated solution evolved in time. In particular, the Lax
spectrum was computed using the Floquet-Fourier-Hill method [70]. Finally, to
reduce the number of simulations needed a double-averaging technique was used
to compute the kurtosis. That is, at each time we compute the spatial average
of $|q|^{2n}$, $n=1,2$, and then we compute the ensemble average.

Figure A7: Same as Fig. A6, but for $m=0.7$.

Figure A8: Same as Fig. A6, but for $m=0.9$.

### Further simulation results

For consistency, the numerical experiments presented in this work related to
the kurtosis have $\epsilon=1/20$, and $\sigma=10^{-2}$ fixed for all values
of $m\in(0,1)$ considered. Table 1 provides key parameters used in the
simulations as well as some of the key quantities related to the kurtosis
experiments. Moreover, Fig. A5 illustrates the kurtosis for $m=0.05$, $0.3$,
$0.5$, $0.7$, $0.8$, and $0.9$ as a function of time. The gray region is the
envelope of 100 ensemble realizations. The red curves are ensemble averages of
$M$ realizations. Notice for each $m$ the ensemble average of the kurtosis
settles at large times (see Table 1 for approximate thermalization times).

To further illustrate thermalization of the SG at large times several plots of
the numerically computed solution of fNLS in the semiclassical limit are
provided. In each case, the left column depicts the time evolution of the dn
potential (A4) while the right column depicts the time evolution of the dn
potential (A4) with a small complex-valued random perturbation. Further, the
top row depicts a density plot while the bottom row depicts the solution at a
particular time. In Fig A6 we have $m=0.05$ which is a small sinusoidal
perturbation of the constant background studied in [34]. In Fig. A8 $m=0.7$,
and in Fig. A8 $m=0.9$. Notice the numerically computed solutions with noise
are spatially homogeneous at large times.

## References

  * [1] M. J. Ablowitz and H. Segur, J. Fluid Mech. 92 ,691–715 (1979)
  * [2] G. P. Agrawal, Nonlinear fiber optics (Academic Press, New York, 2007)
  * [3] E. Infeld and G. Rowlands, Nonlinear waves, solitons and chaos (Cambridge University Press, 2000)
  * [4] L. P. Pitaevskii and S. Stringari, Bose-Einstein condensation (Clarendon Press, Oxford, 2003)
  * [5] M. J. Ablowitz and H. Segur, Solitons and the inverse scattering transform (SIAM, Philadelphia, 1981)
  * [6] S. P. Novikov, S. V. Manakov, L. P. Pitaevskii, and V. E. Zakharov, Theory of solitons: The inverse scattering transform (Plenum, 1984)
  * [7] L D Faddeev and L A Takhtajan, Hamiltonian methods in the theory of solitons (Springer, Berlin, 1987)
  * [8] E. D. Belokolos, A. I. Bobenko, V. Z. Enol’skii, A. R. Its and V. B. Matveev, Algebro-geometric approach to nonlinear integrable equations (Springer-Verlag, Berlin, 1994)
  * [9] F. Gesztesy and H. Holden, Soliton equations and their algebro–geometric solutions vol. I (Cambridge University Press, 2003)
  * [10] M. J. Ablowitz, B. Prinari, and A. D. Trubatch, Discrete and continuous nonlinear Schrödinger systems (Cambridge University Press, 2004)
  * [11] C. S. Gardner, J. M Greene, M. D. Kruskal, and R. M. Miura, Phys. Rev. Lett., 19, 1095 (1967)
  * [12] P. D. Lax, Commun. Pure and Appl. Math., 21, 467–490 (1968)
  * [13] V. E. Zakharov and A. B. Shabat, Sov. Phys. JETP 34, 62–69 (1972)
  * [14] V.E. Zakharov and A. B. Shabat, Sov. Phys. JETP 37, 823 (1973)
  * [15] A. R. Its and V. P. Kotlyarov, Dokl. Akad. Nauk Ukr. 10, 965–968 (1976). English translation: arxiv 1401.4445 [nlin.si]
  * [16] Y.-C. Ma and M. J. Ablowitz, Stud. Appl. Math. 63, 113–158 (1981)
  * [17] M. G. Forest and J. E. Lee, pp. 35–69 (Springer, New York, 1986)
  * [18] D. W. McLaughlin and E. A. Overman II, Surv. Appl. Math. Vol. I, Ed. J. P. Keller, D. W. McLaughlin, G. P. Papanicolaou (Plenum, 1995)
  * [19] G. Biondini and G. Kovačič, J. Math. Phys. 55, 031506 (2014)
  * [20] V. E. Zakharov, Stud. Appl. Math., 122, 219–234 (2009)
  * [21] T. Congy, G. A. El, G. Roberti, A. Tovbis, S. Randoux, and P. Suret, Phys. Rev. Lett. 132, 207201 (2024)
  * [22] G. A. El, Phys. Lett. A 311, 374–383 (2003)
  * [23] G. A. El, J. Stat. Mech., 114001 (2021)
  * [24] P. Suret, S. Randoux, A. Gelash, D. Agafontsev, B. Doyon, and G. El Phys. Rev. E 109, 061001 (2024)
  * [25] A. Costa, A. R. Osborne, D. T. Resio, S. Alessio, E. Chriví, E. Saggese, K. Bellomo, and C. E. Long, Phys. Rev. Lett., 113, 108501 (2014)
  * [26] I. Redor, E. Barthélemy, H. Michallet, M. Onorato, and N. Mordant, Phys. Rev. Lett., 122, 214502 (2019)
  * [27] P. Suret, A. Tikan, F. Bonnefoy, F. Copie, G. Ducrozet, A. Gelash, G. Prabhudesai, G. Michel, A. Cazaubiel, E. Falcon, G. A. El, and S. Randoux, Phys. Rev. Lett., 125, 264101 (2020)
  * [28] T. Leduque, E. Barthélemy, H. Michallet, J. Sommeria, N. Mordant arXiv:2405.14733 [physics.flu-dyn]
  * [29] G. Marcucci, D. Pierangeli, A. J. Agranat, R-K Lee, E. DelRe, and C. Conti, Nat. Commun., 10, 5090 (2019)
  * [30] A. Schwache and F. Mitschke, Phys. Rev. E 55, 7720–7725 (1997).
  * [31] F. Mitschke, G. Steinmeyer, and A. Schwache, Physica D 96, 251–258 (1996),
  * [32] I. Bouchoule and J. Dubail, J. Stat. Mech.2022, 014003 (2022).
  * [33] A. Gelash, D. Agafontsev, V. E. Zakharov, G. El, S. Randoux, and P. Suret, Phys. Rev. Lett. 123, 234102 (2019)
  * [34] A. Gelash and D. S. Agafontsev, Phys. Rev. E 98, 042210 (2018)
  * [35] V. E. Zakharov, Soviet Physics JETP, 33, 538–541 (1971)
  * [36] G. A. El and A. M. Kamchatnov, Phys. Rev. Lett., 95, 204101 (2005)
  * [37] G. A. El and A. Tovbis, Phys. Rev. E 101, 052207 (2020)
  * [38] G. Biondini, X.-D. Luo, J. Oregero and A. Tovbis, Adv. Math. 429, 109188 (2023)
  * [39] G Roberti, G El, A Tovbis, F Copie, P Suret, S Randoux Phys. Rev. E 103, 042205 (2021)
  * [40] D. S. Agafonstev and V. E. Zakharov, Nonlinearity, 28, 2791 (2015)
  * [41] V. E. Zakharov and D.S. Agafontsev, Nonlinearity 29, 3551 (2016)
  * [42] G Roberti, G El, A Tovbis, F Copie, P Suret, S Randoux, Physical Review E 103, 042205 (2021)
  * [43] D. Agafontsev, T.Congy, S. Randoux, P. Suret, G. Roberti and G. El to be submitted
  * [44] The Lax spectrum $\Sigma({\mathcal{L}})$ is defined as the set of values $z\in\mathbb{C}$ such that the focusing ZS scattering problem admits solutions ${\mathbf{v}}(x,z)$ that are bounded for all $x\in\mathbb{R}$.
  * [45] F. W. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University Press, 2010)
  * [46] P. F. Byrd and M. D. Friedman, Handbook of Elliptic Integrals for Engineers and Scientists, (Springer-Verlag, 1971)
  * [47] J. Satsuma and N Yajima, Progr. Theor. Phys. Suppl. 55 284 (1974)
  * [48] Since $\Sigma({\mathcal{L}})$ is invariant when the potential $q(x,t)$ evolves in time according to the fNLS Eq. (1), we will often omit the time dependence from all quantities below.
  * [49] F. Gesztesy and R. Weikard, Acta Math. 181, 63–108 (1998)
  * [50] Please see the appendix for further details.
  * [51] G. Biondini and J. Oregero, Stud. Appl. Math. 145, 325 (2020)
  * [52] G. Biondini, J. Oregero, and A. Tovbis, J. Spectr. Theory 12, 939–992 (2022)
  * [53] A. Tovbis and F. Wang, J. Phys. A: Math. Theor., 55, 424006 (2022)
  * [54] Note that $\int\nolimits_{0}^{K_{m}}\mathop{\rm dn}\nolimits(x;m)\,\mathrm{d}x=\pi/2$, which ensures that $\int_{0}^{1}\varphi(z)\,\mathrm{d}|z|=1$.
  * [55] Note that $f(\eta)$ is a “reduced” DOS. The “full” DOS for the bound state potential is given by $\tilde{f}(z)=f(\eta)\delta(\xi)$, where $\xi=\mathop{\rm Re}\nolimits z$ and $\delta(x)$ is the Dirac delta-function, see [21].
  * [56] M. Onorato, S. Residori, U. Bortolozzo, A. Montina and F.T. Arecchi, Physics Reports 528, 47-89 (2013)
  * [57] P. Walczak, S. Randoux, and P. Suret, Phys. Rev. Lett., 114, 143903 (2015)
  * [58] F. Copie, S. Randoux and P. Suret, Reviews in Physics 5, 100037 (2020)
  * [59] D. Agafontsev, A. Gelash, S. Randoux and P. Suret, Stud. Appl. Math. 152, 810-834 (2024)
  * [60] S. Venakides, Trans. AMS, 301, 189–226 (1987)
  * [61] S. Kamvissis, K. D. T-R McLaughlin, and P. Miller, Semiclassical soliton ensembles for the focusing nonlinear Schrödinger equation (Princeton University Press, 2003)
  * [62] N. M Ercolani, S. Jin, C. D. Levermore, and W. D. MacEvoy, Jr., Int. Math. Res. Not., 2003, 2529–2564 (2003)
  * [63] A. Tovbis, S. Venakides, and X. Zhou, Commun. Pure and Appl. Math., 57, 877–985 (2004)
  * [64] E. G. Shurgalina and E. N. Pelinovsky, Phys. Lett. A, 380 2049–2053 (2016)
  * [65] M. Girotti, T. Grava, R. Jenkins, K. D. T-R McLaughlin, and A. Minakov, Commun. Pure and Appl. Math., 76, 3233–3299 (2023)
  * [66] W. Magnus and S. Winkler, Hill’s equation, Dover, New York, 1966
  * [67] S. Trillo, G. Deng, G. Biondini, M. Klein, G. F. Clauss, A. Chabchoub, and M. Onorato, Phys. Rev. Lett., 117, 144102
  * [68] A. R. Osborne and L. Bergamasco, Phys. D, 18, 26–46 (1986)
  * [69] H. Yoshida, Phys. Lett. A 150, 262–268 (1990)
  * [70] B. Deconinck and J. N. Kutz, J. Comput. Phys. 219, 296–321 (2006)


11institutetext:  Computer Science and Artificial Intelligence Laboratory,
MIT, USA 22institutetext: Centre for Medical Image Computing, University
College London, UK 33institutetext: Biomedical Imaging Group Rotterdam,
Erasmus MC, Netherlands 44institutetext: Laboratory of NeuroImaging,
University of Southern California, USA 55institutetext: Center for Imaging of
Neurodegenerative Diseases, UCSF, USA 66institutetext: Department of Radiology
and Nuclear Medicine, VU Medical Centre, Netherlands
77institutetext: Dementia Research Centre, UCL Institute of Neurology, UK
77email<EMAIL_ADDRESS>
# TADPOLE Challenge: Accurate Alzheimer’s disease prediction through
crowdsourced forecasting of future data

Răzvan V. Marinescu 1122 Neil P. Oxtoby 22 Alexandra L. Young 22 Esther E.
Bron 33 Arthur W. Toga 44 Michael W. Weiner 55 Frederik Barkhof 3366 Nick C.
Fox 77 Polina Golland 11 Stefan Klein 33 Daniel C. Alexander 22

###### Abstract

The Alzheimer’s Disease Prediction Of Longitudinal Evolution (TADPOLE)
Challenge compares the performance of algorithms at predicting the future
evolution of individuals at risk of Alzheimer’s disease. TADPOLE Challenge
participants train their models and algorithms on historical data from the
Alzheimer’s Disease Neuroimaging Initiative (ADNI) study. Participants are
then required to make forecasts of three key outcomes for ADNI-3 rollover
participants: clinical diagnosis, Alzheimer’s Disease Assessment Scale
Cognitive Subdomain (ADAS-Cog 13), and total volume of the ventricles – which
are then compared with future measurements. Strong points of the challenge are
that the test data did not exist at the time of forecasting (it was acquired
afterwards), and that it focuses on the challenging problem of cohort
selection for clinical trials by identifying fast progressors. The submission
phase of TADPOLE was open until 15 November 2017; since then data has been
acquired until April 2019 from 219 subjects with 223 clinical visits and 150
Magnetic Resonance Imaging (MRI) scans, which was used for the evaluation of
the participants’ predictions. Thirty-three teams participated with a total of
92 submissions. No single submission was best at predicting all three
outcomes. For diagnosis prediction, the best forecast (team _Frog_), which was
based on gradient boosting, obtained a multiclass area under the receiver-
operating curve (MAUC) of 0.931, while for ventricle prediction the best
forecast (team _EMC1_), which was based on disease progression modelling and
spline regression, obtained mean absolute error of 0.41% of total intracranial
volume (ICV). For ADAS-Cog 13, no forecast was considerably better than the
benchmark mixed effects model (_BenchmarkME_), provided to participants before
the submission deadline. Further analysis can help understand which input
features and algorithms are most suitable for Alzheimer’s disease prediction
and for aiding patient stratification in clinical trials. The submission
system remains open via the website: https://tadpole.grand-challenge.org/

###### Keywords:

Alzheimer’s Disease, Future prediction, Community Challenge

##  1 Introduction

Accurate prediction of the onset of Alzheimer’s disease (AD) and its
longitudinal progression is important for care planning and for patient
selection in clinical trials. Early detection will be critical in the
successful administration of disease modifying treatments during
presymptomatic phases of the disease prior to widespread brain damage, i.e.
when pathological amyloid and tau accumulate [1]. Moreover, accurate
prediction of the evolution of subjects at risk of Alzheimer’s disease will
help to select homogeneous patient groups for clinical trials, thus reducing
variability in outcome measures that can obscure positive effects on subgroups
of patients who were at the right stage to benefit.

Figure 1: TADPOLE Challenge design. Participants are required to train a
predictive model on a training dataset (D1 and/or others) and make forecasts
for different datasets (D2, D3) by the submission deadline. Evaluation will be
performed on a test dataset (D4) that is acquired after the submission
deadline.

Several approaches for predicting AD-related target variables (e.g. clinical
diagnosis, cognitive/imaging biomarkers) have been proposed which leverage
multimodal biomarker data available in AD. Traditional longitudinal approaches
based on statistical regression model the relationship of the target variables
with other known variables, such as clinical diagnosis [2], cognitive test
scores [3], or time to conversion between diagnoses [4]. Another approach
involves supervised machine learning techniques such as support vector
machines, random forests, and artificial neural networks, which use pattern
recognition to learn the relationship between the values of a set of
predictors (biomarkers) and their labels (diagnoses). These approaches have
been used to discriminate AD patients from cognitively normal individuals [5],
and for discriminating at-risk individuals who convert to AD in a certain time
frame from those who do not [6]. The emerging approach of disease progression
modelling [7, 8] aims to reconstruct biomarker trajectories or other disease
signatures across the disease progression timeline, without relying on
clinical diagnoses or estimates of time to symptom onset. Such models show
promise for predicting AD biomarker progression at group and individual
levels. However, previous evaluations within individual publications are not
systematic and reliable because: (1) they use different data sets or subsets
of the same dataset, different processing pipelines and different evaluation
metrics and (2) over-training can occur due to heavy use of popular training
datasets. Currently we lack a comprehensive comparison of the capabilities of
these methods on standardised tasks relevant to real-world applications.

Community challenges have consistently proven effective in moving forward the
state-of-the-art in technology to address specific data-analysis problems by
providing platforms for unbiased comparative evaluation and incentives to
maximise performance on key tasks. For Alzheimer’s disease prediction in
particular, previous challenges include the CADDementia challenge [9] which
aimed to identify clinical diagnosis from MRI scans. A similar challenge, the
“International challenge for automated prediction of MCI from MRI data“ [10]
asked participants to predict diagnosis and conversion status from extracted
MRI features of subjects from the ADNI study [11]. Yet another challenge, The
Alzheimer’s Disease Big Data DREAM Challenge [12], asked participants to
predict cognitive decline from genetic and MRI data. However, most of these
challenges have not evaluated the ability of algorithms to predict clinical
diagnosis and other biomarkers at future timepoints and largely used training
data from a limited set of modalities. The one challenge that asked
participants to estimate a biomarker at future timepoints (cognitive decline
in one of the DREAM sub-challenges) used only genetic and cognitive data for
training, and aimed to find genetic loci that could predict cognitive decline.
Therefore, standardised evaluation of algorithms needs to be done on biomarker
prediction at future timepoints, with the aim of improving clinical trials
through enhanced patient stratification.

The Alzheimer’s Disease Prediction Of Longitudinal Evolution (TADPOLE)
Challenge aims to identify the data, features and approaches that are most
predictive of future progression of subjects at risk of AD. The challenge
focuses on forecasting the evolution of three key AD-related domains: clinical
diagnosis, cognitive decline, and neurodegeneration (brain atrophy). In
contrast to previous challenges, our challenge is designed to inform clinical
trials through identification of patients most likely to benefit from an
effective treatment, i.e., those at early stages of disease who are likely to
progress over the short-to-medium term (defined as 1-5 years). Since the test
data did not exist at the time of forecast submissions, the challenge provides
a performance comparison substantially less susceptible to many forms of
potential bias than previous studies and challenges. The design choices were
published [13] before the test set was acquired and analysed. TADPOLE also
goes beyond previous challenges by drawing on a vast set of multimodal
measurements from ADNI which support prediction of AD progression.

This article presents the design of the TADPOLE Challenge and outlines
preliminary results.

##  2 Competition Design

The aim of TADPOLE is to predict future outcome measurements of subjects at-
risk of AD, enrolled in the ADNI study. A history of informative measurements
from ADNI (imaging, psychology, demographics, genetics, etc.) from each
individual is available to inform forecasts. TADPOLE participants were
required to predict future measurements from these individuals and submit
their predictions before a given submission deadline. Evaluation of these
forecasts occurred post-deadline, after the measurements had been acquired. A
diagram of the TADPOLE flow is shown in Fig 1.

TADPOLE challenge participants were required to make month-by-month forecasts
of three key biomarkers: (1) clinical diagnosis which is either cognitively
normal (CN), mild cognitive impairment (MCI) or probable Alzheimer’s disease
(AD); (2) Alzheimer’s Disease Assessment Scale Cognitive Subdomain (ADAS-Cog
13) score; and (3) ventricle volume (divided by intra-cranial volume). TADPOLE
forecasts are required to be probabilistic and some evaluation metrics will
account for forecast probabilities provided by participants.

##  3 ADNI data aggregation and processing

TADPOLE Challenge organisers provided participants with a standard ADNI-
derived dataset (available via the Laboratory Of NeuroImaging data archive at
adni.loni.usc.edu) to train algorithms, removing the need for participants to
pre-process the ADNI data or merge different spreadsheets. Software code used
to generate the standard datasets is openly available on
Github111https://github.com/noxtoby/TADPOLE. The challenge data includes: (1)
CSF markers of amyloid-beta and tau deposition; (2) various imaging modalities
such as magnetic resonance imaging (MRI), positron emission tomography (PET)
using several tracers: FDG (hypometabolism), AV45 (amyloid), AV1451 (tau) as
well as diffusion tensor imaging (DTI); (3) cognitive assessments such as
ADAS-Cog 13 acquired in the presence of a clinical expert; (4) genetic
information such as alipoprotein E4 (APOE4) status extracted from DNA samples;
and (5) general demographic information such as age and gender. Extracted
features from this data were merged into a final spreadsheet and made
available online.

The imaging data was pre-processed with standard ADNI pipelines. For MRI
scans, this included correction for gradient non-linearity, B1 non-uniformity
correction and peak sharpening222see http://adni.loni.usc.edu/methods/mri-
analysis/mri-pre-processing. Meaningful regional features such as volume and
cortical thickness were extracted using Freesurfer. Each PET image (FDG, AV45,
AV1451), which consists of a series of dynamic frames, had its frames co-
registered, averaged across the dynamic range, standardised with respect to
the orientation and voxel size, and smoothed to produce a uniform resolution
of 8mm full-width/half-max (FWHM)333see http://adni.loni.usc.edu/methods/pet-
analysis/pre-processing. Standardised uptake value ratio (SUVR) measures for
relevant regions-of-interest were extracted after registering the PET images
to corresponding MR images using SPM5. DTI scans were corrected for head
motion and eddy-current distortion, skull-stripped, EPI-corrected, and finally
aligned to the T1 scans. Diffusion tensor summary measures were estimated
based on the Eve white-matter atlas.

###  3.1 TADPOLE Datasets

In order to evaluate the effect of different methodological choices, we
prepared four “standard“ data sets: the D1 standard training set contains
longitudinal data from the entire ADNI history; the D2 longitudinal prediction
set contains all available data from the ADNI rollover individuals, for whom
challenge participants are asked to provide forecasts; the D3 cross-sectional
prediction set contains a single (most recent) time point and a limited set of
variables from each rollover individual – this represents the information
typically available in a clinical trial; the D4 test set contains visits from
ADNI rollover subjects after 1 Jan 2018, which contain at least one of the
following: diagnostic status, ADAS score, or ventricle volume from MRI – this
dataset did not exist at the time of submitting forecasts. Full demographics
for D1–D4 are given in Table 1.

Measure | D1 | D2 | D3 | D4
---|---|---|---|---
Subjects | 1667 | 896 | 896 | 219
| Cognitively Normal
Number (% total) | 508 (30%) | 369 (41%) | 299 (33%) | 94 (42%)
Visits per subject | 8.3 $\pm$ 4.5 | 8.5 $\pm$ 4.9 | 1.0 $\pm$ 0.0 | 1.0 $\pm$ 0.2
Age | 74.3 $\pm$ 5.8 | 73.6 $\pm$ 5.7 | 72.3 $\pm$ 6.2 | 78.4 $\pm$ 7.0
Gender (% male) | 48% | 47% | 43% | 47%
MMSE | 29.1 $\pm$ 1.1 | 29.0 $\pm$ 1.2 | 28.9 $\pm$ 1.4 | 29.1 $\pm$ 1.1
Converters (% total CN) | 18 (3.5%) | 9 (2.4%) | - | -
| Mild Cognitive Impairment
Number (% total) | 841 (50.4%) | 458 (51.1%) | 269 (30.0%) | 90 (41.1%)
Visits per subject | 8.2 $\pm$ 3.7 | 9.1 $\pm$ 3.6 | 1.0 $\pm$ 0.0 | 1.1 $\pm$ 0.3
Age | 73.0 $\pm$ 7.5 | 71.6 $\pm$ 7.2 | 71.9 $\pm$ 7.1 | 79.4 $\pm$ 7.0
Gender (% male) | 59.3% | 56.3% | 58.0% | 64.4%
MMSE | 27.6 $\pm$ 1.8 | 28.0 $\pm$ 1.7 | 27.6 $\pm$ 2.2 | 28.1 $\pm$ 2.1
Converters (% total MCI) | 117 (13.9%) | 37 (8.1%) | - | 9 (10.0%)
| Alzheimer’s Disease
Number (% total) | 318 (19.1%) | 69 (7.7%) | 136 (15.2%) | 29 (13.2%)
Visits per subject | 4.9 $\pm$ 1.6 | 5.2 $\pm$ 2.6 | 1.0 $\pm$ 0.0 | 1.1 $\pm$ 0.3
Age | 74.8 $\pm$ 7.7 | 75.1 $\pm$ 8.4 | 72.8 $\pm$ 7.1 | 82.2 $\pm$ 7.6
Gender (% male) | 55.3% | 68.1% | 55.9% | 51.7%
MMSE | 23.3 $\pm$ 2.0 | 23.1 $\pm$ 2.0 | 20.5 $\pm$ 5.9 | 19.4 $\pm$ 7.2
Converters (% total AD) | - | - | - | 9 (31.0%)
Table 1: Summary of TADPOLE datasets D1–D4. Each subject has been allocated to
either Cognitively Normal, MCI or AD group based on diagnosis at the first
available visit within each dataset.

##  4 Submissions and evaluation

The challenge had a total of 33 participating teams, who submitted a total of
58 forecasts from D2, 34 forecasts from D3, and 6 forecasts from custom
prediction sets. Table 2 summarises the top-3 winner methods in terms of input
features used, handling of missing data and predictive models: _Frog_ used a
gradient boosting method, which combined many weak predictors to build a
strong predictor; _EMC1_ derived a “disease state“ variable aggregating
multiple features together and then used an SVM and 2D splines for prediction,
while VikingAI used a latent-time parametric model with subject- and feature-
specific parameters – see [14] for full method details. We also describe three
benchmark models which were provided to participants at the start of the
challenge: (i) _BenchmarkLastVisit_ uses the measurement at the last available
visit, (ii) _BenchmarkME-APOE_ uses a mixed effects model with APOE status as
covariate and (iii) BenchmarkSVM uses an out-of-the-box support vector machine
(SVM) and regressor for forecast.

Submission |  Extra† Features |  Nr. of features |  Missing data imputation |  Diagnosis prediction |  ADAS/Vent. prediction
---|---|---|---|---|---
Frog |  most features |  70+420* |  none |  gradient boosting |  gradient boosting
EMC1-Std |  MRI, ASL, cognitive |  250 |  nearest neighbour |  DPM SVM 2D-spline |  DPM 2D-spline
VikingAI-Sigmoid |  MRI, cognitive, tau |  10 |  none |  DPM + ordered logit |  DPM
BenchmarkLastVisit |  - |  3 |  none |  constant model |  constant model
BenchmarkME-APOE |  APOE |  4 |  none |  Gaussian model |  linear mixed effects model
BenchmarkSVM |  age, APOE |  6 |  mean of previous values |  SVM |  support vector regressor
Table 2: Summary of benchmarks and top-3 methods used in the TADPOLE
submissions. DPM – disease progression model. (†) Aside from the three target
biomarkers (*) Augmented features: e.g. min/max, trends, moments.

For evaluation of clinical status predictions, we used similar metrics to
those that proved effective in the CADDementia challenge [9]: (i) the
multiclass area under the receiver operating curve (MAUC); and (ii) the
overall balanced classification accuracy (BCA). For ADAS and ventricle volume,
we used three metrics: (i) mean absolute error (MAE), (ii) weighted error
score (WES) and (iii) coverage probability accuracy (CPA). BCA and MAE focus
purely on prediction accuracy ignoring confidence, MAUC and WES include
confidence, while CPA provides an assessment of the confidence interval only.
Complete formulations for these can be found in Table 3, with detailed
explanations in the TADPOLE design paper [13]. To compute an overall rank, we
first calculated the sum of ranks from MAUC, ADAS MAE and Ventricle MAE for
each submission, and the overall ranking was derived from these sums of ranks.

Formula |  Definitions
---|---
$mAUC=\frac{2}{L(L-1)}\sum_{i=2}^{L}\sum_{j=1}^{i}\hat{A}(c_{i},c_{j})$ |  $n_{i}$, $n_{j}$ – number of points from class $i$ and $j$. $S_{ij}$ – the sum of the ranks of the class $i$ test points, after ranking all the class $i$ and $j$ data points in increasing likelihood of belonging to class $i$, $L$ – number of data points
$BCA=\frac{1}{2L}\sum_{i=1}^{L}\left[\frac{TP}{TP+FN}+\frac{TN}{TN+FP}\right]$ |  $TP_{i}$, $FP_{i}$, $TN_{i}$, $FN_{i}$ – the number of true positives, false positives, true negatives and false negatives for class $i$ $L$ – number of data points
$MAE=\frac{1}{N}\sum_{i=1}^{N}\left|{\tilde{M}_{i}-M_{i}}\right|$ |  $M_{i}$ is the actual value in individual $i$ in future data. $\tilde{M}_{i}$ is the participant’s best guess at $M_{i}$ and $N$ is the number of data points
$WES=\frac{\sum_{i=1}^{N}\tilde{C}_{i}\left|\tilde{M}_{i}-M_{i}\right|}{\sum_{i=1}^{N}\tilde{C}_{i}}$ |  $M_{i}$, $\tilde{M}_{i}$ and $N$ defined as above. $\tilde{C}_{i}=(C_{+}-C_{-})^{-1}$, where $\left[C_{-},C_{+}\right]$ is the 50% confidence interval
$CPA=|ACP-0.5|$ |  actual coverage probability (ACP) - the proportion of measurements that fall within the 50% confidence interval.
Table 3: TADPOLE performance metric formulas and definitions for the terms.

##  5 Results

While full results can be found on the TADPOLE website [14], here we only
include the top-3 winners. Table 4 compiles all metrics for top-3 TADPOLE
forecasts from the D2 prediction set. The best overall performance was
obtained by team _Frog_ , with a clinical diagnosis MAUC of 0.931, ADAS MAE of
4.85 and Ventricle MAE of 0.45. Among the benchmark methods, _BenchmarkME-
APOE_ had the best overall rank of 18, obtaining an MAUC of 0.82, ADAS MAE of
4.75 and Ventricle MAE of 0.57. In terms of diagnosis predictions, _Frog_ had
an overall MAUC score of 0.931. For ADAS prediction, _BenchmarkME-APOE_ had
the best MAE of 4.75. For Ventricle prediction, _EMC1-Std_ had the best MAE of
0.41 and WES of 0.29. In terms of the most accurate confidence interval
estimates, _VikingAI_ achieved the best CPA scores of 0.02 for ADAS and 0.2
for Ventricles.

| Overall | Diagnosis | ADAS | Ventricles (% ICV)
---|---|---|---|---
Submission | Rank | MAUC | BCA | MAE | WES | CPA | MAE | WES | CPA
Frog | 1 | 0.931 | 0.849 | 4.85 | 4.74 | 0.44 | 0.45 | 0.33 | 0.47
EMC1-Std | 2 | 0.898 | 0.811 | 6.05 | 5.40 | 0.45 | 0.41 | 0.29 | 0.43
VikingAI-Sigmoid | 3 | 0.875 | 0.760 | 5.20 | 5.11 | 0.02 | 0.45 | 0.35 | 0.20
BenchmarkME-APOE | 18 | 0.822 | 0.749 | 4.75 | 4.75 | 0.36 | 0.57 | 0.57 | 0.40
BenchmarkSVM | 34 | 0.836 | 0.764 | 6.82 | 6.82 | 0.42 | 0.86 | 0.84 | 0.50
BenchmarkLastVisit | 40 | 0.774 | 0.792 | 7.05 | 7.05 | 0.45 | 0.63 | 0.61 | 0.47
Table 4: Ranked forecasting scores for benchmark models and top-3 TADPOLE
submissions.

##  6 Discussion

In the current work we have outlined the design and key results of TADPOLE
Challenge, which aims to identify algorithms and features that can best
predict the evolution of Alzheimer’s disease. Despite the small number of
converters in the training set, the methods were able to accurately forecast
the clinical diagnosis and ventricle volume, although they found it harder to
forecast cognitive test scores. Compared to the benchmark models, the best
submissions had considerably smaller errors that represented only a small
fraction of the errors obtained by benchmark models (0.42 for clinical
diagnosis MAUC and 0.71 for ventricle volume MAE). For clinical diagnosis,
this suggests that more than half of the subjects originally misdiagnosed by
the best benchmark model (_BenchmarkSVM_) are now correctly diagnosed with the
new methods. Moreover, the results suggest that we do not have a clear winner
on all categories. While team Frog had the best overall submission with the
lowest sum of ranks, for each performance metric individually we had different
winners.

Additional work currently in progress [14] suggests that consensus methods
based on averaging predictions from all participants perform better than any
single individual method. This demonstrates the power of TADPOLE in achieving
state-of-the-art prediction accuracy through crowd-sourcing prediction models.

The TADPOLE Challenge and its preliminary results presented here are of
importance for the design of future clinical trials and more generally may be
applicable to a clinical setting. The best algorithms identified here could be
used for subject selection or stratification in clinical trials, e.g. by
enriching trial inclusion with fast progressors to increase the statistical
power to detect treatment changes. Alternatively, a stratification could be
implemented based on predicted “fast progressors“ and “slow progressors“ to
reduce imbalances between arms. In order to make these models applicable to
clinical settings, application in a clinical sample should be tested outside
ADNI and further validation in a subject population with post-mortem
confirmation would be desirable, as clinical diagnosis of probable AD only has
moderate agreement with gold-standard neuropathological post-mortem diagnosis
(70.9% – 87.3% sensitivity and 44.3% – 70.8% specificity, according to [15]).
We hope such a validation will be possible in the future, with the advent of
neuropathological confirmation in large, longitudinal, multimodal datasets
such as ADNI.

In future work, we plan to analyse which features and methods were most useful
for predicting AD progression, and assess if the results are sufficient to
improve stratification for AD clinical trials. We also plan to evaluate the
impact and interest of the first phase of TADPOLE within the community, to
guide decisions on whether to organise further submission and evaluation
phases.

##  7 Acknowledgements

TADPOLE Challenge has been organised by the European Progression Of
Neurological Disease (EuroPOND) consortium, in collaboration with the ADNI. We
thank all the participants and advisors, in particular Clifford R. Jack Jr.
from Mayo Clinic, Rochester, United States and Bruno M. Jedynak from Portland
State University, Portland, United States for useful input and feedback.

The organisers are extremely grateful to The Alzheimer’s Association, The
Alzheimer’s Society and Alzheimer’s Research UK for sponsoring the challenge
by providing the prize fund and providing invaluable advice into its
construction and organisation. Similarly, we thank the ADNI leadership and
members of our advisory board and other members of the EuroPOND consortium for
their valuable advice and support.

RVM was supported by the EPSRC Centre For Doctoral Training in Medical Imaging
with grant EP/L016478/1 and by the Neuroimaging Analysis Center with grant NIH
NIBIB NAC P41EB015902. NPO, FB, SK, and DCA are supported by EuroPOND, which
is an EU Horizon 2020 project. ALY was supported by an EPSRC Doctoral Prize
fellowship and by EPSRC grant EP/J020990/01. PG was supported by NIH grant
NIBIB NAC P41EB015902 and by grant NINDS R01NS086905. DCA was supported by
EPSRC grants J020990, M006093 and M020533. The UCL-affiliated researchers
received support from the NIHR UCLH Biomedical Research Centre. Data
collection and sharing for this project was funded by the Alzheimer’s Disease
Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01
AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012).
FB was supported by the NIHR UCLH Biomedical Research Centre and the AMYPAD
project, which has received support from the EU-EFPIA Innovative Medicines
Initiatives 2 Joint Undertaking (AMYPAD project, grant 115952). This project
has received funding from the European Union Horizon 2020 research and
innovation programme under grant agreement No 666992.

## References

  * [1] Mehta, D., Jackson, R., Paul, G., Shi, J., Sabbagh, M.:  Why do trials for Alzheimer’s disease drugs keep failing? A discontinued drug perspective for 2010–2015.  Expert opinion on investigational drugs 26(6) (2017) 735
  * [2] Scahill, R.I., Schott, J.M., Stevens, J.M., Rossor, M.N., Fox, N.C.:  Mapping the evolution of regional atrophy in Alzheimer’s disease: unbiased analysis of fluid-registered serial MRI.  Proceedings of the National Academy of Sciences 99(7) (2002) 4703–4707
  * [3] Yang, E., Farnum, M., Lobanov, V., Schultz, T., Raghavan, N., Samtani, M.N., Novak, G., Narayan, V., DiBernardo, A.:  Quantifying the pathophysiological timeline of Alzheimer’s disease.  Journal of Alzheimer’s Disease 26(4) (2011) 745–753
  * [4] Guerrero, R., Schmidt-Richberg, A., Ledig, C., Tong, T., Wolz, R., Rueckert, D., ADNI, et al.:  Instantiated mixed effects modeling of Alzheimer’s disease markers.  NeuroImage 142 (2016) 113–125
  * [5] Klöppel, S., Stonnington, C.M., Chu, C., Draganski, B., Scahill, R.I., Rohrer, J.D., Fox, N.C., Jack Jr, C.R., Ashburner, J., Frackowiak, R.S.:  Automatic classification of MR scans in Alzheimer’s disease.  Brain 131(3) (2008) 681–689
  * [6] Young, J., Modat, M., Cardoso, M.J., Mendelson, A., Cash, D., Ourselin, S., ADNI, et al.:  Accurate multimodal probabilistic prediction of conversion to Alzheimer’s disease in patients with mild cognitive impairment.  NeuroImage: Clinical 2 (2013) 735–745
  * [7] Young, A.L., Oxtoby, N.P., Daga, P., Cash, D.M., Fox, N.C., Ourselin, S., Schott, J.M., Alexander, D.C.:  A data-driven model of biomarker changes in sporadic Alzheimer’s disease.  Brain 137(9) (2014) 2564–2577
  * [8] Lorenzi, M., Filippone, M., Frisoni, G.B., Alexander, D.C., Ourselin, S., Initiative, A.D.N., et al.:  Probabilistic disease progression modeling to characterize diagnostic uncertainty: application to staging and prediction in Alzheimer’s disease.  NeuroImage 190 (2017) 56–68
  * [9] Bron, E.E., Smits, M., Van Der Flier, W.M., Vrenken, H., Barkhof, F., Scheltens, P., Papma, J.M., Steketee, R.M., Orellana, C.M., Meijboom, R., et al.:  Standardized evaluation of algorithms for computer-aided diagnosis of dementia based on structural MRI: the CADDementia challenge.  NeuroImage 111 (2015) 562–579
  * [10] Sarica, A., Cerasa, A., Quattrone, A., Calhoun, V.:  Editorial on special issue: Machine learning on MCI.  Journal of neuroscience methods 302 (2018) 1
  * [11] Weiner, M.W., Veitch, D.P., Aisen, P.S., Beckett, L.A., Cairns, N.J., Green, R.C., Harvey, D., Jack Jr, C.R., Jagust, W., Morris, J.C., et al.:  Recent publications from the Alzheimer’s Disease Neuroimaging Initiative: Reviewing progress toward improved AD clinical trials.  Alzheimer’s & Dementia 13(4) (2017) e1–e85
  * [12] Allen, G.I., Amoroso, N., Anghel, C., Balagurusamy, V., Bare, C.J., Beaton, D., Bellotti, R., Bennett, D.A., Boehme, K.L., Boutros, P.C., et al.:  Crowdsourced estimation of cognitive decline and resilience in Alzheimer’s disease.  Alzheimer’s & dementia: the journal of the Alzheimer’s Association 12(6) (2016) 645–653
  * [13] Marinescu, R.V., Oxtoby, N.P., Young, A.L., Bron, E.E., Toga, A.W., Weiner, M.W., Barkhof, F., Fox, N.C., Klein, S., Alexander, D.C., et al.:  Tadpole challenge: Prediction of longitudinal evolution in Alzheimer’s disease.  arXiv preprint arXiv:1805.03909 (2018)
  * [14] https://tadpole.grand-challenge.org/Results/
  * [15] Beach, T.G., Monsell, S.E., Phillips, L.E., Kukull, W.:  Accuracy of the clinical diagnosis of alzheimer disease at national institute on aging alzheimer disease centers, 2005–2010.  Journal of neuropathology and experimental neurology 71(4) (2012) 266–273