# Self-organized dynamics of a viscous drop with interfacial nematic activity

Mohammadhossein Firouznia  Department of Mechanical and Aerospace Engineering,
University of California San Diego,
9500 Gilman Drive, La Jolla, California, 92093, USA  Center for Computational
Biology, Flatiron Institute, New York, New York, 10010, USA
David Saintillan<EMAIL_ADDRESS>Department of Mechanical and Aerospace
Engineering, University of California San Diego,
9500 Gilman Drive, La Jolla, California, 92093, USA

###### Abstract

We study emergent dynamics in a viscous drop subject to interfacial nematic
activity. Using hydrodynamic simulations, we show how the interplay of
nematodynamics, activity-driven flows and surface deformations gives rise to a
sequence of self-organized behaviors of increasing complexity, from periodic
braiding motions of topological defects to chaotic defect dynamics and active
turbulence, along with spontaneous shape changes and translation. Our findings
recapitulate qualitative features of experiments and shed light on the
mechanisms underpinning morphological dynamics in active interfaces.

Living materials are characterized by their ability to continuously transform
chemical energy into mechanical work at the microscale. When confined to
deformable surfaces, these active materials demonstrate a myriad of dynamic
behaviors. Biological surfaces also often exhibit intrinsic degrees of freedom
that correspond to in-plane order (such as nematic or polar), which facilitate
long-range hydrodynamic interactions, resulting in the emergence of self-
organized spatiotemporal patterns. This is a fundamental feature of various
biological systems, including the cell cortex and confluent eukaryotic cells,
and plays a crucial role in their functional properties. Here, we focus on
systems with nematic symmetry. Experimental observations indicate the
emergence of nematic order during cytokinesis [1, 2]. Nematic alignment has
also been evidenced in different stages of tissue morphogenesis when
individual cells exhibit a preferred elongation axis [3, 4, 5].

Biological processes from subcellular to multicellular scales often occur on
self-deforming surfaces with various mechanical properties. Examples include
shape changes during polarization, migration and division in cells, and apical
constrictions in epithelial morphogenesis [6, 7]. Understanding these
morphological dynamics necessitates theoretical models that account for the
three-way coupling between surface deformations, in-plane order, and flow [8,
9, 10]. Recent studies have probed the mechanisms driving shape changes in
cells and epithelial tissues, advancing our understanding of these processes
[11, 12, 13].

The central motif of these biological systems has also been utilized to create
synthetic soft active materials [14, 15, 16]. In pioneering work, Keber et al.
[14] assembled a shape-shifting lipid vesicle by encapsulating a film of
microtubules and kinesin motors to its inner surface. Understanding the
interplay between orientational order, activity-induced flow, and substrate
geometry has been the subject of several studies thereafter. Various models
have probed the role of surface curvature on the dynamics of topological
defects in active fluids confined to rigid surfaces of various topologies [17,
18, 19, 20, 21], yet the role of interfacial deformations and their coupling
to bulk flow remains poorly understood.

In this Letter, we report on the spontaneous dynamics of a viscous drop driven
out of equilibrium due to interfacial nematic activity. We show that the
interplay between the flow inside and outside the drop, surface transport of
the nematic field and surface deformations, gives rise to a sequence of self-
organized behaviors and symmetry-breaking phenomena of increasing complexity.
Our results recapitulate the qualitative features of experiments [14], both in
the small and finite deformation regimes. Under small deformations, the
dynamics is characterized by the braiding motion of topological defects around
the drop, giving rise to well-known braiding patterns at different activity
levels. The asymmetry induced under finite deformations results in
translational motion of the drop. Under strong activity, a transition to
active turbulence is also observed.

We model a viscous drop occupying volume $V^{-}$ and suspended in another
viscous fluid $V^{+}$ [Fig. 1(a)]. The interface $\partial V$ is defined by a
smooth surface $\bm{r}=\bm{r}(s^{1},s^{2})\in\mathbb{R}^{3}$ parameterized by
coordinates $(s^{1},s^{2})$. Tangent vectors $\bm{g}_{i}=\partial_{i}\bm{r}$
($\partial_{i}\\!\coloneqq\partial/\partial s^{i}$, $i=1,2$), along with the
unit normal
$\smash{\bm{v}={(\bm{g}_{1}\times\bm{g}_{2})}/{|\bm{g}_{1}\times\bm{g}_{2}|}}$,
form a local coordinate system with surface metric tensor
$g_{ij}=\bm{g}_{i}\bm{\cdot}\bm{g}_{j}$. The drop is initially spherical with
radius $R$. A monolayer of active nematic particles is constrained to
$\partial V$ and drives the system out of equilibrium by inducing flows inside
and outside the drop and by causing deformations.

Figure 1: (a) Schematic of the model system. (b) Snapshots of the director
field and scalar order parameter shown as a colormap during: (left) braiding
motion under small deformations, (middle) braiding motion with shape-shifting
behavior, and (right) chaotic regime with finite deformations. Also see Videos
1–3 of the SM [22]. (c) Tetrahedral (left) and planar (right) defect
arrangements. (d,e) Evolution of the average angle $\alpha$ for two systems
with $\mathrm{Pe}=0.3$ (d) and $\mathrm{Pe}=3$ (e). The red and blue lines
denote planar ($120^{\circ}$) and tetrahedral ($109.5^{\circ}$)
configurations, respectively. $(\mathrm{Ca},\,\ell_{c})=(0.01,\,0.32)$ in
(d,e).

We use a coarse-grained representation of the nematic field in terms of the
tensor $\bm{Q}=Q^{ij}\,\bm{g}_{i}\bm{g}_{j}$ where
$Q^{ij}=S(n^{i}n^{j}-g^{ij}/2)$. Here, $\bm{n}=n^{i}\bm{g}_{i}$
($n^{i}n_{i}=1$) denotes the local nematic director, and $S\in[0,1]$ is the
scalar order parameter characterizing the strength of alignment. The $\bm{Q}$
tensor evolves by the nematodynamic equation [23, 21, 9]

$\displaystyle\mathrm{D}_{t}Q^{ij}=\dfrac{1}{\Gamma}H^{ij}+\zeta\widetilde{U}^{ij},\qquad\bm{x}\in\partial
V.$ (1)

$\mathrm{D}_{t}Q^{ij}$ is a co-rotational material derivative embodying
advection and rotation of the nematic by the surface flow with velocity
$\bm{u}=u^{i}\bm{g}_{i}+u_{n}\bm{v}$:

$\begin{split}\mathrm{D}_{t}Q^{ij}=&~{}\partial_{t}Q^{ij}+u^{k}\,\nabla_{k}Q^{ij}+u_{n}\,(C^{j}_{k}Q^{ik}+C^{i}_{k}Q^{kj})\\\
&\qquad+\omega_{n}\,(\epsilon^{ik}Q^{j}_{k}+\epsilon^{jk}Q^{i}_{k}),\end{split}$
(2)

where $\nabla_{k}$ is the covariant derivative with respect to $s^{k}$,
$C_{ij}=-\partial_{j}\partial_{i}\bm{x}\bm{\cdot}\bm{n}$ is the curvature
tensor, $\omega_{n}=\frac{1}{2}\epsilon^{ij}\nabla_{i}u_{j}$ is the normal
vorticity, and $\epsilon^{ij}=\bm{v}\bm{\cdot}(\bm{g}^{i}\times\bm{g}^{j})$ is
the Levi-Civita tensor. The molecular tensor $H^{ij}=-\delta\mathcal{F}/\delta
Q^{ij}$ describes orientational relaxation in the nematic monolayer, with
$\Gamma$ the rotational viscosity. It derives from the Landau–de Gennes free
energy,

$\mathcal{F}\\!=\\!\\!\int_{\partial
V}\\!\\!\\!\\!\mathrm{d}A\Bigl{(}\frac{k_{s}}{2}\\!\left[d\,Q_{ij}Q^{ij}\\!+\\!c(Q_{ij}Q^{ij})^{2}\right]\\!+\\!\frac{k_{e}}{2}\nabla_{i}Q_{jk}\nabla^{i}Q^{jk}\\!\Bigr{)},$
(3)

which accounts for short-range and elastic interactions in the nematic
monolayer [24, 25]. Here, $k_{s}$ and $k_{e}$ are steric and elastic
phenomenological constants, and $d=a+\frac{b}{3}S+\frac{c}{6}S^{2}$ where $a$,
$b$ and $c$ are normalized thermotropic parameters. We assume that the
characteristic size of a nematic particle is small compared to the local
radius of curvature of the interface, and thus neglect any coupling between
the extrinsic curvature and the nematic tensor in Eq. (3). Finally, alignment
by the flow is captured by $\smash{\zeta\widetilde{U}^{ij}}$ where $\zeta$ is
the flow alignment parameter or Bretherton constant [26], and
$\widetilde{U}^{ij}$ is the traceless strain-rate tensor:

$\widetilde{U}^{ij}={U}^{ij}-\tfrac{1}{2}U^{k}_{k}g^{ij},\,\,\,{U}^{ij}=\tfrac{1}{2}(\nabla^{i}u^{j}+\nabla^{j}u^{i})+C^{ij}u_{n}.$
(4)

Neglecting inertial effects and gravity, the flow inside and outside the drop
is governed by the Stokes equations:

${\mu^{\pm}}\nabla^{2}\bm{u}^{\pm}-\bm{\nabla}p^{\pm}=\mathbf{0},\quad\bm{\nabla}\bm{\cdot}\bm{u}^{\pm}=0,\quad~{}\bm{x}\in
V^{\pm}.$ (5)

The velocity is continuous across $\partial V$ and vanishes far away. The
nematic particles exert an active surface stress $\bm{T}^{a}=\xi\,\bm{Q}$ on
their surrounding, leading to fluid motion and deformations. The constant
$\xi$ captures the biochemical activity, with $\xi<0$ for extensile systems as
in the experiments of Keber et al. [14]. The local force balance along the
tangential and normal directions on the drop surface $\partial V$ reads

$\displaystyle f^{h,j}+\xi\,\nabla_{i}Q^{ij}=0,$ (6)
$\displaystyle{f^{h}_{n}-\gamma\,C_{k}^{k}-\,\xi\,C_{ij}\,Q^{ij}=0,}$ (7)

where $\gamma$ is the uniform surface tension. The jump in hydrodynamic
tractions across $\partial V$ is
$\bm{f}^{h}=\bm{v}\bm{\cdot}[\bm{T}^{h}\big{]}_{-}^{+}$, where
$\bm{T}^{h}=-p\bm{I}+\mu\big{(}\bm{\nabla}\bm{u}+\bm{\nabla}\bm{u}^{\mathrm{T}}\big{)}$
is the Newtonian stress tensor. We solve Eqs. (1)–(7) numerically using a
custom spectral boundary integral solver [27, 28], where all physical
variables such as shape, velocity, and $\bm{Q}$ tensor are represented as
truncated series of spherical surface harmonics (see Supplemental Material
(SM) [22] for details).

Dimensional analysis of the governing equations yields four dimensionless
groups in addition to $\zeta$. The active capillary number
$\mathrm{Ca}={|\xi|\,R}/{\gamma}$ compares the strength of active stresses to
surface tension and governs the magnitude of deviations from the spherical
shape [Fig. 1(b)]. The active Péclet number
$\mathrm{Pe}={|\xi|\,\Gamma}/{(k_{s}\mu^{+}R^{2})}$ characterizes the strength
of convective vs relaxational fluxes in Eq. (1). The balance between short-
range and elastic interactions in the nematic monolayer defines a
dimensionless coherence length $\ell_{c}=R^{-1}\sqrt{k_{e}/k_{s}}$, or
effective distance over which topological defects affect the nematic field on
the scale of the drop. Finally, the viscosity ratio between the inner and
outer fluids is $\lambda=\mu^{-}/\mu^{+}$. We found that varying $\lambda$ has
little effect on the dynamics, and all results shown here are for $\lambda=1$.

Figure 2: (a) Spatiotemporal trajectories $(X,Y,t)$ of projected defects for a
drop with $(\mathrm{Pe},\ell_{c})=(0.3,0.16)$ during golden braiding pattern
(see SM for details [22]). (b) Diagrams of golden (triangle) and silver
(square) braids. (c) Defect trajectories in 3D for active drops with
$\ell_{c}=0.16$ in the golden (left), silver (middle), and mix (right)
braiding regimes at different values of $\mathrm{Pe}$. The time interval
between the first and last data points is $\Delta t=50/f_{1}$, where $f_{1}$
corresponds to the peak frequency in the FFT spectrum of $\alpha$. (d) FTBE of
defect trajectories, and (e) average active power $\bar{P}_{a}$ for a
spherical active drop, as functions of $\mathrm{Pe}$ and $\ell_{c}$. The
markers identify golden, silver, and mix braiding regimes. $\mathrm{Ca}=0.01$
in all cases.

First, we focus on the regime of $\mathrm{Ca}\ll 1$ to isolate the effect of
deformations. In this regime, the behavior of the system is governed by two
parameters: the activity level captured by $\mathrm{Pe}$, and nematic
elasticity captured by $\ell_{c}$. Active stresses drive the system out of
equilibrium while nematic elasticity tends to stabilizes it. Indeed, as
nematic elasticity becomes stronger, the coherence length increases
($\ell^{2}_{c}\propto k_{e}$), resulting in a more constrained arrangement of
the topological defects that repel one another. Given the spherical topology
of the interface, the net topological charge of the nematic field is fixed at
$+2$, and for small to moderate activity levels ($0<\mathrm{Pe}\lesssim 5$),
the system exhibits four $+1/2$ defects [Fig. 1(b), left].

Under small activity levels, the defects remain stationary and occupy the
vertices of a regular tetrahedron (Fig. S2 [22]). Beyond a critical value
$\mathrm{Pe}_{c}$, a transition occurs to an unsteady regime marked by the
periodic motion of the defects (see Video 1 [22]). Each defect is set into
motion by the activity-induced flow, which influences the nematic field. In
agreement with experimental observations [14], the defect arrangement
oscillates between tetrahedral and planar configurations, as shown in Fig.
1(c-e), where $\alpha=1/6\sum_{i<j}^{4}\alpha_{ij}$ denotes the average of all
pairwise angles $\alpha_{ij}$ between the four defects.

The spatiotemporal defect trajectories can be analyzed and interpreted using
the concept of braids [29]. We first map the defect trajectories from the
spherical surface of the drop to its two-dimensional mid-plane using the
stereographic projection, where one reference defect is taken to be at the
sphere pole and therefore mapped to infinity [30, 31]. When projected along a
specific direction, the remaining three defect trajectories display a sequence
of crossings that defines a braid $\bm{b}$ [Fig. 2(a)]. Each crossing event is
expressed in terms of elementary generators
$\sigma_{i},~{}i\in\\{1,\dots,n-1\\},$ of the $n$-particle braid group with
$n=3$ [29]. The generator $\sigma_{i}$ denotes the clockwise exchange of
defect $i$ with defect $i+1$ on the projection line, while $\sigma_{i}^{-1}$
corresponds to their counterclockwise exchange. The complexity of a given
braid and its mixing efficiency can be quantified by the exponential
stretching rate of material curves in the associated flow. This is determined
through the braid’s topological entropy and finite-time braiding exponent
(FTBE) for periodic and aperiodic trajectories, respectively [32, 33, 34] (see
SM for details [22]).

For activity levels close to $\mathrm{Pe}_{c}$ within the unsteady regime, the
defect motion is periodic and characterized by a single time scale, as
evidenced by the presence of a single peak in the fast Fourier transform (FFT)
of the average angle $\alpha$ (Fig. S4 [22]). The defect trajectories follow a
specific pattern known as the golden braid [Fig. 2(a,b)], described by the
braid $\bm{B}_{1}={(\bm{b}_{1})}^{k}$, where
$\bm{b}_{1}=\sigma_{1}\,\sigma_{2}^{-1}\,\sigma_{1}\,\sigma_{1}\,\sigma_{2}\,\sigma_{1}^{-1}\,\sigma_{2}^{-1}\,\sigma_{2}^{-1}$.
Its topological entropy is given by $h(\bm{b}_{1})\\!=6\log{\phi_{1}}$, where
$\phi_{1}=(1+\sqrt{5})/2$ denotes the golden ratio [35]. As $\mathrm{Pe}$ is
further increased, nonlinear effects become more pronounced, leading to the
emergence of a second time scale. This is manifested in the FFT spectrum of
$\alpha$, which exhibits two distinct peaks (Fig. S5 [22]). Eventually, the
braiding pattern undergoes a transition from the golden braid to the silver
braid $\bm{B}_{2}={(\bm{b}_{2})}^{k}$, where
$\bm{b}_{2}=\sigma_{1}\,\sigma_{2}\,\sigma_{1}^{-1}\,\sigma_{2}^{-1}\,\sigma_{1}\,\sigma_{2}^{-1}$.
Its topological entropy is $h(\bm{b}_{2})=2\log{\phi_{2}}$, where
$\phi_{2}=1+\sqrt{2}$ is the silver ratio [36]. Schematic diagrams of the
golden and silver braids are shown in Fig. 2(b). The efficiency of a periodic
braid in increasing entropy can be quantified by its topological entropy by
generator (TEPG) [36]. Comparing the TEPG of the golden and silver braids
shows that $\mathrm{TEPG}(\bm{b}_{1})\approx 1.23\mathrm{TEPG}(\bm{b}_{2})$,
i.e., the transition to the silver braid with increasing $\mathrm{Pe}$ causes
a decrease in TEPG.

As the activity level is further increased, a transition from periodic to
aperiodic dynamics is observed. This regime, referred to as mix braiding, is
characterized by the absence of a specific braiding pattern, with defect
trajectories displaying irregular and chaotic behavior. The FFT spectrum of
$\alpha$ in this regime is broad and displays multiple peaks, a signature of
highly nonlinear dynamics (Fig. S6 [22]). Figure 2(c) shows 3D defect
trajectories in drops with golden, silver, and mix braiding patterns,
providing visual evidence of the increasing complexity with higher activity.
To characterize the different braiding patterns and quantify their
complexities in the non-periodic case, we calculate the FTBE of defect
trajectories as a function of $\mathrm{Pe}$ and $\ell_{c}$ in Fig. 2(d) [34,
37]. The FTBE generally shows an increasing trend with respect to
$\mathrm{Pe}$, indicating that higher activity levels lead to more complex
dynamics. However, the transition from golden to silver braiding is
accompanied by a decrease in FTBE, consistent with the decrease in TEPG
discussed above. This reduction in the FTBE is particularly intriguing as it
occurs despite the heightened activity.

To elucidate this effect, we analyze the drop’s energetics in the unsteady
regime. The active power expended by the surface nematic is defined as
$\smash{{P}_{a}(t)=\int_{{}_{\partial
V}}\bm{f}^{a}\bm{\cdot}\bm{u}\,\mathrm{d}s}$, where
$\smash{\bm{f}^{a}=\xi\nabla_{i}Q^{ij}\bm{g}_{j}-\,\xi C_{ij}\,Q^{ij}\bm{v}}$
is the interfacial active traction and $\bm{u}$ the interfacial velocity. As
shown in Fig. 2(e), the time-averaged active power, denoted by
$\bar{P}_{a}=\langle P_{a}\rangle_{t}$, is an ascending function of activity
except during the transition from the golden to silver braid, where it
suddenly drops. This indicates that the increase in TEPG with $\mathrm{Pe}$ in
the golden braiding regime comes at a higher energetic cost, and suggests that
the spontaneous transition to the silver braid is energetically favorable.
Note that, under small deformations, the active power is dissipated primarily
through viscous effects within the bulk, with negligible contributions from
capillarity (see SM for details [22]).

We highlight the stabilizing effect of nematic elasticity as captured by
$\ell_{c}$. The critical Péclet number $\mathrm{Pe}_{c}$ for the transition
from equilibrium to periodic braiding rises from $\mathrm{Pe}_{c}\approx 0.2$
to $1.2$ as the coherence length is varied from $\ell_{c}=0.16$ to $0.32$. The
onset of aperiodic defect motions (mix braiding regime) is also delayed under
larger coherence lengths. According to Fig. 2(d,e), for a given $\mathrm{Pe}$,
both the FTBE and the active power are consistently lower at the larger
coherence length. Indeed, the leading effect of nematic elasticity is to repel
topological defects, which fosters more organized and constrained defect
dynamics and delays the transition to chaotic motion.

Allowing for finite deformations ($\mathrm{Ca}>0$) further increases the
complexity of the dynamics. In this regime, the active drop undergoes
spontaneous shape changes in the form of breathing motions as the defects
traverse its surface [Fig. 1(b), middle] (See Video 2 [22]). These
deformations, in turn, impact the nematodynamics, establishing a three-way
feedback loop between shape, nematic field, and flow. Notably, we observe an
increase in $\mathrm{Pe}_{c}$ at higher capillary numbers, indicating that
elevated activity levels are required for defects to overcome the energy
barriers induced by deformations.

For the smaller coherence length $\ell_{c}=0.16$, we observe that under
moderate values of $\mathrm{Pe}>\mathrm{Pe}_{c}$, four defects exhibit
braiding motion similar to the behavior observed in the small deformation
limit. Concurrently, the drop undergoes spontaneous shape changes under the
influence of active stresses. At long times, we find that the drop eventually
reaches an equilibrium state with a steady deformed shape and defect
configuration. This behavior is only observed up to a second critical Péclet
number $\mathrm{Pe}_{eq}>\mathrm{Pe}_{c}$, beyond which the system transitions
to a chaotic regime characterized by the rapid creation and annihilation of
defects around the drop [Fig. 1(b), right] (also see Video 3 and Fig. S9
[22]). The newly formed defects emerge as pairs with $\pm 1/2$ topological
charges, maintaining a constant net topological charge. The dynamics in this
regime resembles the active turbulence previously observed in nematic
materials on flat and curved surfaces [38, 39, 40, 21, 41]. For the larger
coherence length $\ell_{c}=0.32$, the system reaches equilibrium over long
times for all explored Péclet numbers ($0<\mathrm{Pe}<12$) due to the
stabilizing effect of nematic elasticity.

Figure 3: Time-averaged translational velocity as a function of $\mathrm{Pe}$
for $\mathrm{Ca}=0.5$. Blue and brown lines correspond to $\ell_{c}=0.16$ and
$0.32$, respectively. Filled markers indicate cases where more than four
defects were observed. Insets show snapshots of the drop corresponding to the
boxed markers.

Nonlinearities introduced by finite deformations amplify the asymmetry in the
nematic field and defect configuration. This asymmetry gives rise to a net
translational component in the velocity field. The volume-averaged velocity of
the drop of volume $V_{d}$ is calculated as
$\smash{\bm{U}(t)=\frac{1}{V_{d}}\int_{{V^{-}}}\\!\bm{u}\,\mathrm{d}\upsilon}$,
and we denote by $\bar{U}=|\langle\bm{U}\rangle_{t}|$ the magnitude of its
time average. We observe that the translational velocity is orders of
magnitude larger compared to the regime of $\mathrm{Ca}\ll 1$, underscoring
the role of deformations in breaking the system’s symmetry. Figure 3 shows
$\bar{U}$ vs $\mathrm{Pe}$ for $\mathrm{Ca}=0.5$. For $\ell_{c}=0.16$,
$\bar{U}$ first exhibits a sharp increase with $\mathrm{Pe}$ when four defects
are present. It then reaches a plateau as $\mathrm{Pe}$ is further increased
and the drops enters the chaotic regime where additional defects are created.
Similar trends are observed for $\ell_{c}=0.32$, although the increase in
$\bar{U}$ with $\mathrm{Pe}$ is delayed due to the stabilizing effect of
nematic elasticity, which promotes more symmetric shapes.

Using numerical simulations, we have analyzed the dynamics of a deformable
viscous drop subject to interfacial nematic activity. Our results highlight
the complex interplay of nematodynamics, active fluid flows, and interfacial
mechanics, and point to a wide range of emergent dynamics depending on the
importance of active stresses relative to viscous and capillary stresses. In
the low deformation regime, the dynamics is characterized by four +1/2 defects
whose trajectories can be described in terms of braids of increasing
complexity. At finite capillary numbers, breathing deformations and
translational motion emerge, as well as an active turbulent regime with
continuous generation and annihilation of defect pairs. Our observations are
all consistent with past experiments on active vesicles [14]. We note that the
analysis of braiding motions only provides qualitative information on the
system’s mixing efficiency, and is limited to regimes with four topological
defects, where treating defects as material points is a reasonable
approximation. A more in-depth characterization of the organizing role of
activity-induced flows may rely on an analysis of Lagrangian coherent
structures and associated finite-time Lyapunov exponents, which have recently
been applied to identify flow attractors and repellers in bulk active nematics
[42, 43] as well as during embryogenesis [44]. Finally, we note that the
system studied here may serve as a simplified model for active living systems,
such as cells or organoids. For that purpose, various model extensions may be
desirable, including accounting for the role of elastic stresses, of extrinsic
curvature coupling, or of chemical cues, which are instrumental in regulating
the cell cortex [45, 46, 47].

###### Acknowledgements.

The authors thank J.-L. Thiffeault for insightful conversations on braids and
topological entropy, and S. H. Bryngelson for his contributions to the
computational framework [28]. This work was partially funded by National
Science Foundation Grants CBET-1934199 and DMS-2153520.

## References

  * Spira  _et al._ [2017] F. Spira, S. Cuylen-Haering, S. Mehta, M. Samwer, A. Reversat, A. Verma, R. Oldenbourg, M. Sixt, and D. W. Gerlich, Cytokinesis in vertebrate cells initiates by contraction of an equatorial actomyosin network composed of randomly oriented filaments, eLife 6, e30867 (2017).
  * Reymann  _et al._ [2016] A.-C. Reymann, F. Staniscia, A. Erzberger, G. Salbreux, and S. W. Grill, Cortical flow aligns actin filaments to form a furrow, eLife 5, e17807 (2016).
  * Guillamat  _et al._ [2022] P. Guillamat, C. Blanch-Mercader, G. Pernollet, K. Kruse, and A. Roux, Integer topological defects organize stresses driving tissue morphogenesis, Nat. Mater. 21, 588 (2022).
  * Saw  _et al._ [2017] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux, Topological defects in epithelia govern cell death and extrusion, Nature 544, 212 (2017).
  * Duclos  _et al._ [2017] G. Duclos, C. Erlenkämper, J.-F. Joanny, and P. Silberzan, Topological defects in confined populations of spindle-shaped cells, Nat. Phys. 13, 58 (2017).
  * Salbreux  _et al._ [2012] G. Salbreux, G. Charras, and E. Paluch, Actin cortex mechanics and cellular morphogenesis, Trends Cell Biol. 22, 536 (2012).
  * Martin and Goldstein [2014] A. C. Martin and B. Goldstein, Apical constriction: themes and variations on a cellular mechanism driving morphogenesis, Development 141, 1987 (2014).
  * Mietke  _et al._ [2019] A. Mietke, F. Jülicher, and I. F. Sbalzarini, Self-organized shape dynamics of active surfaces, Proc. Natl. Acad. Sci. U.S.A. 116, 29 (2019).
  * Salbreux  _et al._ [2022] G. Salbreux, F. Jülicher, J. Prost, and A. Callan-Jones, Theory of nematic and polar active fluid surfaces, Phys. Rev. Res. 4, 033158 (2022).
  * Al-Izzi and Morris [2023] S. C. Al-Izzi and R. G. Morris, Morphodynamics of active nematic fluid surfaces, J. Fluid Mech. 957, A4 (2023).
  * Metselaar  _et al._ [2019] L. Metselaar, J. M. Yeomans, and A. Doostmohammadi, Topology and morphology of self-deforming active shells, Phys. Rev. Lett. 123, 208001 (2019).
  * Ruske and Yeomans [2021] L. J. Ruske and J. M. Yeomans, Morphology of active deformable 3D droplets, Phys. Rev. X 11, 021001 (2021).
  * Khoromskaia and Salbreux [2023] D. Khoromskaia and G. Salbreux, Active morphogenesis of patterned epithelial shells, eLife 12, e75878 (2023).
  * Keber  _et al._ [2014] F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, and A. R. Bausch, Topology and dynamics of active nematic vesicles, Science 345, 1135 (2014).
  * Guillamat  _et al._ [2018] P. Guillamat, Ž. Kos, J. Hardoüin, J. Ignés-Mullol, M. Ravnik, and F. Sagués, Active nematic emulsions, Sci. Adv. 4, eaao1470 (2018).
  * Weirich  _et al._ [2019] K. L. Weirich, K. Dasbiswas, T. A. Witten, S. Vaikuntanathan, and M. L. Gardel, Self-organizing motors divide active liquid droplets, Proc. Natl. Acad. Sci. U.S.A. 116, 11125 (2019).
  * Zhang  _et al._ [2016] R. Zhang, Y. Zhou, M. Rahimi, and J. J. De Pablo, Dynamic structure of active nematic shells, Nat. Commun. 7, 13483 (2016).
  * Shankar  _et al._ [2017] S. Shankar, M. J. Bowick, and M. C. Marchetti, Topological sound and flocking on curved surfaces, Phys. Rev. X 7, 031039 (2017).
  * Green  _et al._ [2017] R. Green, J. Toner, and V. Vitelli, Geometry of thresholdless active flow in nematic microfluidics, Phys. Rev. Fluids 2, 104201 (2017).
  * Nitschke  _et al._ [2018] I. Nitschke, M. Nestler, S. Praetorius, H. Löwen, and A. Voigt, Nematic liquid crystals on curved surfaces: a thin film limit, Proc. Math. Phys. Eng. Sci. 474, 20170686 (2018).
  * Pearce  _et al._ [2019] D. Pearce, P. W. Ellis, A. Fernandez-Nieves, and L. Giomi, Geometrical control of active turbulence in curved topographies, Phys. Rev. Lett. 122, 168002 (2019).
  * Sup [2022] (2022), see Supplemental Material available at (link to be inserted by publisher) for additional details on the model, simulation method and analysis, and for videos of simulations corresponding to Fig. 1(b).
  * Giomi [2015] L. Giomi, Geometry and topology of turbulence in active nematics, Phys. Rev. X 5, 031003 (2015).
  * Nitschke  _et al._ [2020] I. Nitschke, S. Reuther, and A. Voigt, Liquid crystals on deformable surfaces, Proc. Math. Phys. Eng. Sci. 476, 20200313 (2020).
  * Kralj  _et al._ [2011] S. Kralj, R. Rosso, and E. G. Virga, Curvature control of valence on nematic shells, Soft Matter 7, 670 (2011).
  * Bretherton [1962] F. P. Bretherton, The motion of rigid particles in a shear flow at low Reynolds number, J. Fluid Mech. 14, 284 (1962).
  * Firouznia  _et al._ [2023] M. Firouznia, S. H. Bryngelson, and D. Saintillan, A spectral boundary integral method for simulating electrohydrodynamic flows in viscous drops, J. Comput. Phys. 489, 112248 (2023).
  * Firouznia  _et al._ [2022] M. Firouznia, S. H. Bryngelson, and D. Saintillan, Spectral Boundary Integral Solver for Electrohydrodynamic Flows in Viscous Drops, https://github.com/mfirouzn/EHD_Drop_3D (2022).
  * Thiffeault [2022] J.-L. Thiffeault,  _Braids and Dynamics_ (Springer Nature, 2022).
  * Fadell and Buskirk [1962] E. Fadell and J. V. Buskirk, The braid groups of $E^{2}$ and $S^{2}$, Duke Math. J. 29, 243 (1962).
  * Birman [1974] J. S. Birman,  _Braids, Links, and Mapping Class Groups_ (Princeton University Press, 1974).
  * Boyland  _et al._ [2000] P. L. Boyland, H. Aref, and M. A. Stremler, Topological fluid mechanics of stirring, J. Fluid Mech. 403, 277 (2000).
  * Fathi  _et al._ [1979] A. Fathi, F. Laundenbach, and V. Poénaru, Travaux de Thurston sur les surfaces, Asterisque 66–67, 1 (1979).
  * Budišić and Thiffeault [2015] M. Budišić and J. Thiffeault, Finite-time braiding exponents, Chaos 25 (2015).
  * D’Alessandro  _et al._ [1999] D. D’Alessandro, M. Dahleh, and I. Mezic, Control of mixing in fluid flow: A maximum entropy approach, IEEE Trans. Autom. Control. 44, 1852 (1999).
  * Finn and Thiffeault [2011] M. D. Finn and J.-L. Thiffeault, Topological optimization of rod-stirring devices, SIAM Review 53, 723 (2011).
  * Thiffeault and Budišić [2021] J. Thiffeault and M. Budišić, Braidlab: A software package for braids and loops (2013–2021), version 3.2.5, arXiv:1410.0849 [math.GT] .
  * Guillamat  _et al._ [2017] P. Guillamat, J. Ignés-Mullol, and F. Sagués, Taming active turbulence with patterned soft interfaces, Nat. Commun. 8, 564 (2017).
  * Ellis  _et al._ [2018] P. W. Ellis, D. J. G. Pearce, Y. Chang, G. Goldsztein, L. Giomi, and A. Fernandez-Nieves, Curvature-induced defect unbinding and dynamics in active nematic toroids, Nat. Phys. 14, 85 (2018).
  * Lin  _et al._ [2021] S. Lin, W. Zhang, D. Bi, B. Li, and X. Feng, Energetics of mesoscale cell turbulence in two-dimensional monolayers, Commun. Phys. 4, 21 (2021).
  * Doostmohammadi  _et al._ [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Active nematics, Nat. Commun. 9, 3246 (2018).
  * Serra  _et al._ [2023] M. Serra, L. Lemma, L. Giomi, Z. Dogic, and L. Mahadevan, Defect-mediated dynamics of coherent structures in active nematics, Nat. Phys. 7, 1366 (2023).
  * Sinigaglia  _et al._ [2023] C. Sinigaglia, F. Braghin, and M. Serra, Optimal control of short-time attractors in active fluid flows, arXiv preprint arXiv:2305.00193 (2023).
  * Serra  _et al._ [2020] M. Serra, S. Streichan, M. Chuai, C. J. Weijer, and L. Mahadevan, Dynamic morphoskeletons in development, Proc. Natl. Acad. Sci. U.S.A. 117, 11444 (2020).
  * Lessey  _et al._ [2012] E. C. Lessey, C. Guilluy, and K. Burridge, From mechanical force to RhoA activation, Biochem. 51, 7420 (2012).
  * Bement  _et al._ [2015] W. M. Bement, M. Leda, A. M. Moe, A. M. Kita, M. E. Larson, A. E. Golding, C. Pfeuti, K.-C. Su, A. L. Miller, A. B. Goryachev, _et al._ , Activator–inhibitor coupling between Rho signalling and actin assembly makes the cell cortex an excitable medium, Nat. Cell Biol. 17, 1471 (2015).
  * Bischof  _et al._ [2017] J. Bischof, C. A. Brand, K. Somogyi, I. Májer, S. Thome, M. Mori, U. S. Schwarz, and P. Lénárt, A cdk1 gradient guides surface contraction waves in oocytes, Nat. Commun. 8, 849 (2017).


# Heterogeneous Image-based Classification Using Distributional Data Analysis

Alec Reinhardt1111Corresponding author: Alec Reinhardt, Department of
Biostatistics, UT MD Anderson Cancer Center, 1400 Pressler St., Floor 4,
FCT4.6000, Houston, TX 77030. E-mail<EMAIL_ADDRESS>Newsha
Nikzad2 Raven J. Hollis2 Galia Jacobson2 Millicent A Roach2 Mohamed Badawy3
Peter Chul Park4 Laura Beretta5 Prasun K Jalal6 David T. Fuentes7 Eugene J.
Koay2,8 and Suprateek Kundu1
1Department of Biostatistics, MD Anderson Cancer Center
2 Cancer Physics and Engineering Laboratory, MD Anderson Cancer Center
3Department of Pediatrics, UT Southwestern Medical Center
4Department of Radiation Oncology, UC Davis Comprehensive Cancer Center
5Department of Molecular and Cellular Oncology, MD Anderson Cancer Center
6Department of Hepatology, Baylor College of Medicine
7Department of Imaging Physics, MD Anderson Cancer Center
8Department of Radiation Oncology, MD Anderson Cancer Center

(August 28, 2024)

Abstract

Diagnostic imaging has gained prominence as potential biomarkers for early
detection and diagnosis in a diverse array of disorders including cancer.
However, existing methods routinely face challenges arising from various
factors such as image heterogeneity. We develop a novel imaging-based
distributional data analysis (DDA) approach that incorporates the probability
(quantile) distribution of the pixel-level features as covariates. The
proposed approach uses a smoothed quantile distribution (via a suitable basis
representation) as functional predictors in a scalar-on-functional quantile
regression model. Some distinctive features of the proposed approach include
the ability to: (i) account for heterogeneity within the image; (ii)
incorporate granular information spanning the entire distribution; and (iii)
tackle variability in image sizes for unregistered images in cancer
applications. Our primary goal is risk prediction in Hepatocellular carcinoma
that is achieved via predicting the change in tumor grades at post-diagnostic
visits using pre-diagnostic enhancement pattern mapping (EPM) images of the
liver. Along the way, the proposed DDA approach is also used for case versus
control diagnosis and risk stratification objectives. Our analysis reveals
that when coupled with global structural radiomics features derived from the
corresponding T1-MRI scans, the proposed smoothed quantile distributions
derived from EPM images showed considerable improvements in sensitivity and
comparable specificity in contrast to classification based on routinely used
summary measures that do not account for image heterogeneity. Given that there
are limited predictive modeling approaches based on heterogeneous images in
cancer, the proposed method is expected to provide considerable advantages in
image-based early detection and risk prediction.

Keywords: Cancer imaging; distributional data analysis; scalar-on-function
regression; quantile regression

##  1 Introduction

Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer
and a significant global health concern. It is the sixth most common cancer
and the fourth leading cause of cancer-related deaths worldwide. One of the
key factors contributing to its high mortality rate is the lack of early
detection. Early detection plays a pivotal role in improving liver cancer
prognosis and increasing treatment success rates. Patients diagnosed at an
early stage have the opportunity to undergo curative treatments such as
resection, ablative therapies, and liver transplantation, while those
diagnosed at a later stage typically only qualify for palliative systemic
treatments with limited effectiveness. Consequently, the 5-year survival rate
surpasses 70% for individuals with early-stage HCC, whereas it falls below 5%
for those diagnosed at advanced stages. Surveillance for HCC is recommended in
at-risk patients, including those with cirrhosis of all etiologies, and
certain populations with chronic HBV infection. Guidelines recommend HCC
surveillance with abdominal ultrasound (US) with or without serum alpha-
fetoprotein (AFP) measurement [36] every 6 months. Other diagnostic imaging
modalities for surveillance include computed tomography (CT) and magnetic
resonance imaging (MRI), and ultrasound (US) imaging. However, the sensitivity
of the surveillance may potentially vary with the imaging modality [32].

Although medical imaging has made rapid advances and has emerged as one of the
most promising tools for early diagnosis and risk detection, there are
important limitations. These include both MRI-based over surveillance due to
the frequent detection of benign vascular lesions, and delayed HCC diagnosis
due to the under-staging of small or early malignant lesions in approximately
25% of cases [21, 26]. A potential reason for these inadequacies is the
inherent noise in liver MRI scans, which is often unaccounted for. Another
major diagnostic and prognostic challenge of rising concern is the
heterogeneity of pre-cancerous and cancerous lesions. Despite technological
advances, lesion detection, characterization of patterns of enhancement, and
monitoring remain sub-optimal. While the Liver Imaging Reporting and Data
System (LI-RADS) has been developed as a standardized classification system
for liver nodules (depending on lesion size, arterial phase enhancement,
capsular enhancement, washout, and threshold growth to deduce a nodule’s
probability of malignancy, and so on), risk stratification is hindered by
significant heterogeneity within LI-RADS categories, especially for LR-3 and
LR-4 lesions of which 38% and 74% are HCC, respectively [4, 24]. A recent
meta-analysis showed that there was significant heterogeneity in sensitivity
of US-based surveillance (21%-89%) across studies included in the meta-
analysis [32], which can be potentially attributed to the above factors and
highlights a limitation of imaging-based surveillance. Fuethermore, in other
domains such as pancreatic cancer, standardized scoring system such as LI-RADS
may not always be available. This exacerbates difficulties for imaging-based
risk stratification, and the standard of care may resort to biopsy or other
invasive procedures for a definitive diagnostic assessment. Therefore, there
is critical need for new minimally invasive imaging platforms and accompanying
statistical approaches that can render risk stratification with high
sensitivity and specificity based on imaging data collected at earlier disease
stages.

In this article, we propose novel statistical and machine learning models for
early detection and risk prediction in cancer imaging studies, motivated by
HCC application. In particular, our goal is to develop a novel distributional
data analysis (DDA) based approach leveraging quantile distributions [22]
across imaging pixels to be used for classification or prediction. This
article makes several important contributions.

First, we use a novel type of liver image known as enhancement pattern mapping
(EPM) that are derived from T1-MRI scans and have higher contrast-to-noise
ratio. EPM is a novel pixel-based signal analysis technique that quantifies
the difference in enhancement over time of a given pixel in the liver compared
to either a patient-specific or population-based normal liver model, providing
a measurement and visualization of how different the signal is over an entire
volume of interest. Previous literature [26] has shown that EPM derived from
CT images improves the contrast-to-noise (CNR) ratio enhancement for lesion
detection in hepatobiliary malignancy. In contrast to this previous research,
we use EPM features derived from multi-phase MRI scans in order to evaluate
the predictive potential of a novel image-based biomarker. Secondly from a
methodological perspective, we propose a novel distributional data analysis
(DDA) approach that uses the entire quantile distribution of pixel-level
imaging measurements for classification and risk prediction. Such an approach
extends the DDA literature proposed for regression [35] to cancer imaging
involving high-dimensional heterogeneous images. By utilizing full information
embedded in the probability distribution of granular pixel-level measurements,
the proposed approach is able to account for the heterogeneity within the
image that may arise from the differences of histopathological properties of
lesional, peri-lesional and healthy tissues. Therefore, this approach provides
considerable advantages over routinely used approaches based on summary
measures. Another important advantage of the DDA based method is that one does
not require the images to be registered across samples and the methodology can
handle lesions of varying sizes that is often encountered in cancer imaging
applications. Third, we develop a novel scalar-on-functional quantile
classification methodology that leverages penalization techniques to achieve
sparsity and model parsimony in the selected features. Fourth, we focus on
several important analysis goals including case vs control diagnosis and risk
stratification, along with the critical goal of early detection that involves
identifying lesions presenting the highest risk of longitudinal progression
based on data from the initial disease stages. These goals are described in
more detail below.

Our initial objective involves disease classification, aiming to distinguish
between cases and controls using whole liver images with and without lesion
mask information. The subsequent task is risk stratification, focusing on
classifying the grade of HCC liver lesions (mild versus aggressive). Following
this, we shift our attention to the primary objective of early detection,
which holds significant clinical and translational importance by identifying
lesions with high potential for longitudinal progression early on, thus
improving survival outcomes. Specifically, in aim 3, we employ the proposed
scalar-on-functional quantile classification approach to identify lesions
likely to exhibit aggressive longitudinal progression based on imaging data
from the pre-diagnostic visit. We predict the tumor grade at the post-
diagnostic visit using training images from the pre-diagnostic visit and
evaluate performance accordingly. In the above classification analyses, we
compare the performance of different combinations of functional quantile
distributional features and structural radiomics features as well as other
image summaries. Collectively, these combined analyses establish a
comprehensive prediction framework based on pixel-level heterogeneous imaging
data, capable of addressing gaps in literature. In particular, existing
approaches primarily rely on summary imaging features without accounting for
heterogeneity with the image, resulting in pitfalls.

Furthermore, we conduct thorough numerical comparisons demonstrating the
notable benefits of integrating functional quantile distribution features with
structural radiomics features for classification. This surpasses the
effectiveness of using functional quantile feature or structural radiomics
features by themselves or even routinely used summary measures derived from
the image. Our meticulous numerical experiments reveal that even with a
limited to moderate number of functional quantlet features in the
classification model, the proposed method achieves sensitivity levels
exceeding 80% and specificity nearing 95% for early detection. Additionally,
we observe consistent performance across varying lesion sizes and despite
small to moderate training sample sizes. These findings clearly underscore the
potential of the proposed functional quantile distributions extracted from EPM
images in cancer diagnosis and risk prediction, especially when combined with
structural radiomics features. A schematic diagram of the proposed analytical
pipeline is presented in Figure 1.

Figure 1: A schematic of the proposed analytical pipeline. The leftmost panel
captures the 5 different phases of the T1-MRI liver scan that are all
registered to a common template for a given subject. Subsequently the EPM
images are computed from these T1-MRI scans, and then empirical quantile
distributions followed by smoothed quantile distributions (based on quantlet
basis expansions) were computed. The arterial phase of T1-MRI is used for
segmentation of tumor and also used to compute the structural radiomics
features. These two types of features, along with demographics are used for
three classification tasks that include early detection, risk stratification,
and diagnosis (case vs controls).

##  2 Related Literature

There is a growing statistical and machine learning literature for image-based
modeling as summarized below.

Analysis using vectorized pixels: The earliest approaches for image-based
analysis relied on mass univariate regression analysis [9] that may not
account for dependence across imaging pixels. Later, more sophisticated Markov
Random Field (MRF) pixel-wise regression approaches have been proposed that
are able to account for the spatial dependence [29, 19]. While useful, pixel-
level approaches described above often encounter the curse of dimensionality
due to tens of thousands of imaging pixels in a single image. Recently,
tensor-based approaches have been proposed [15, 12] that simultaneously
alleviate the curse of dimensionality and account for the spatial information
in the images. However, the above approaches were often designed for
neuroimaging data where the images are required to be registered across
subjects. While registration of brain images is standard in neuroimaging
literature, established pipelines for image registration across samples are
typically not readily available for cancer imaging studies due to various
factors such as varying size, composition, and other characteristics of the
organs of interest. Therefore, the above sophisticated methods are not
directly generalizable to cancer imaging applications that is our focus.

Analysis based on summary radiomics features: Summary features such as mean
pixel intensity [3] as well as more informative radiomic features are
extracted from high-dimensional radiological images [17] via established
algorithms. These features capture different global features of the image such
as first order and second order statistics, as well as image texture. Radiomic
analyses have been applied within oncology to improve diagnosis and prognosis
[16, 1] with the aim of delivering precision medicine. The radiomic workflow
involves curation of clinical and imaging data and is a step-wise process
involving image pre-processing, tumor segmentation, feature extraction, model
development, and validation [18]. Such a pipeline requires input from
individuals in many disciplines and at multiple levels, resulting in an overly
supervised approach that is potentially susceptible to subject-specific bias.
Moreover, radiomics features collapse information across the whole image or
regions of interest (ROI) that leads to information loss, and they may fail to
capture the heterogeneity within the image, which is a key (but often
overlooked) concern [5].

Scalar-on-function models: Methods for functional data analysis [27] have
become ubiquitous with the growth of recent technologies that are able to
generate high-dimensional functional data. Although the vast majority of
literature has focused on one-dimensional functional curves [23], recent
literature has started investigating models involving more complex types of
functional data such as images. For example, neuroimaging analysis using
entire brain images as covariates (scalar-on-image regression) for prediction
of a continuous outcome have gained in popularity [8, 20]. In these types of
models, the image is treated as a functional object. Typical scalar-on-image
regression approaches account for hundreds of thousands of pixels via low-
dimensional representation of images, for instance principal components or
tensor methods. For example, approaches involving functional principal
component analysis (FPCA) [37] often assume that the components driving
variability in the images are related to the outcome that may not always be
practical. An advantage of these dimension reduction methods is that they do
not require the image to be registered across samples. However, FPCA
approaches only use a subset of principal components resulting in information
loss and potentially compromising predictive accuracy. Alternative approaches
involving basis expansion representations have been proposed for scalar-on-
image regression [34, 28] and provide a desirable avenue to preserve the
spatial properties of the image. Unfortunately, such basis expansion methods
typically require images to be registered across samples that may not feasible
in cancer imaging studies. Moreover, basis expansion approaches typically do
not have in-built mechanisms for low dimensional feature extraction and
therefore may suffer from the curse of dimensionality.

Statistical literature on distributional data analysis: Distributional data
analysis (DDA) can be considered a specific form of functional data analysis,
where the covariate probability (or quantile) distribution is used as the
functional predictor. This type of approach requires repeated measures from a
given observed covariate, which is increasingly common in recent studies with
growing volumes of data. In such studies, each covariate in a data set is
associated with its own probability distribution that represents the
heterogeneity in the features of interest, which is subsequently used in
modeling approaches. Such studies may involve user-defined histograms to be
used as predictors in scalar-on-function regression [2], or compositional data
analysis (CoDA) [7] that refers to another group of methods to model
continuously measured wearable data to analyze samples of densities [30].
Distribution-on-distribution regression models were suggested by [10].
Additional work on DDA involve research by [11] who used subject-specific
quantile functions to capture the distributional nature of wearable data, and
[22] who used a quantile-function representation of accelerometry data to
predict health outcomes. While the use of DDA approaches is growing rapidly,
these techniques have not yet been adopted for mainstream cancer imaging
applications accounting for heterogeneity, to our knowledge.

##  3 Methods and Materials

###  3.1 Description of Liver Imaging Dataset

The study population included patients with cirrhosis presenting to the
Hepatology Clinic between 2012 and 2020 at a single tertiary care hospital who
were prospectively followed. Patients underwent HCC surveillance with the
instructional practice of contrast-enhanced MRI every six months. Patients
with cirrhosis were included in the study if they had at least two consecutive
contrast-enhanced MRIs, utilizing liver protocol for HCC surveillance.
Patients with cirrhosis and HCC diagnosed from baseline imaging were excluded.
The dataset subsequently comprised of  48 individuals with LR-3 and LR-4 liver
lesions (cases), and 97 healthy subjects without any lesions (controls).
Among the case subjects, we excluded certain subjects and/or visits based on
quality control for the EPM images and lesion segmentations. After exclusion,
imaging and clinical data on 38 case subjects were available for analysis, out
of which 35 subjects had imaging data for both timepoints. Among the  73 case
images (35 subject $\times$ 2 visits + 3 subjects $\times$ 1 visit), 91
distinct lesion areas were identified including multiple lesions for a  subset
of subjects. These 91 lesions include all pre-diagnostic lesions and post-
diagnostic lesions but exclude minuscule-sized lesions with about 10 or less
pixels. Longitudinal differences in tumor characteristics and grades were
observed between the pre-diagnostic and post-diagnostic images for 35 subjects
with imaging data for both timepoints. In addition to the MRI scans and
corresponding lesion and liver masks, the LIRADS score for each lesion were
also available that ranged from 2-5. Arterial phase hyperenhancement of the
lesion, venous and delayed phases washout, in addition to the lesion’s size
and growth pattern were the criteria used to assign category codes based on
LI-RADS version 2018 guidelines. Spatially distinct lesioned areas are
identified using a clustering algorithm. When multiple distinct lesioned areas
are identified, we separate them for processing when the nearest points lie at
least 5 mm apart. Figure 2 provides visual illustrations.

####  3.1.1 Enhancement Pattern Mapping (EPM) Images

The original T1w-MRI scans were processed to obtain voxel-wise EPM data that
corresponds to the blood perfusion in the liver tissue. EPM is a three-
dimensional, voxel-based algorithm that was earlier used for quantitative
image analysis of liver lesions based on CT scans [26], which we now adapt for
MRI scans. The generalized enhancement pattern, such as quantitative intensity
changes throughout multi-phase MRI due to uptake and washout of contrast
materials was acquired from the registered multi-phase MRI scans. A normal
liver enhancement curve was developed through fitting the MRI intensity values
over time within user-selected ROIs, sampled uniformly across cirrhotic liver
parenchyma from the same patient. The root-mean-square deviation (RMSD) for
each voxel was computed by averaging the squares of differences between the
generalized normal liver intensity and the voxel intensity across all time
points, and subsequently taking the square root of the average. The RMSD
values of each voxel was mapped to their original MRI coordinates, and these
values comprised the EPM image maps used for our analysis. The normal liver
enhancement curve was obtained by fitting MRI intensity values sampled from
the normal liver ROIs over the period of contrast phases by a piece-wise
smooth function, where each piece was a second-order polynomial. The EPM
algorithm was implemented numerically using MATLAB [14].

###  3.2 Extraction of radiomics features

We engineer and extract novel functional radiomics features from the EPM
images, and combine them with standard structural radiomics features from the
MRI scans in our analysis.

Functional Quantlet Features based on EPM: We computed functional quantlet
features that capture the information embedded in the distribution of EPM
values across all pixels in the lesion and/or the liver (depending on the
analysis goals). These functional features encompass a much richer set of
information spanning the entire distribution compared to summary statistics
typically reported in existing pipelines such as [33]. The functional quantlet
features rely on a basis expansion that enables smoothing and potential
interpolation, thus reducing the impact of unequal-sized lesions. Further,
they enable a parsimonious representation via a near-lossless property that
ensures that only a moderate number of these quantlets are able to capture
almost all of the information in the quantile distribution. The full
mathematical details for deriving the quantlet features are described in the
next section.

Structural radiomics features based on T1-MRI: In addition to EPM images,
structural radiomics features were also extracted based on the arterial phase
of MRI scans. For controls, each ROI (6mm in diameter) was manually defined
over background liver parenchyma, and subsequently the structural radiomics
features were extracted from these ROIs using the Pyradiomics package [33].
For cases, the structural radiomics features were computed separately for each
lesion as well as for an independent non-lesional (pre-defined) ROI. Both
liver masks and lesion masks were used for defining ROIs, where the lesion
segmentation was conducted based on the arterial phase of the T1w-MRI scan.
When computing the structural radiomics features for cases, we did not
differentiate between multiple lesions with the same LIRADS score.

From a total of 108 radiomics features extracted from suitable ROIs, we
included four categories of structural radiomics information: 24 Gray Level
Co-occurrence Matrix (GLCM) features, 16 Gray Level Run Length Matrix (GLRLM)
features, 16 Gray Level Size Zone Matrix (GLSZM) features, and 14 shape-based
features. We excluded variables representing first-order structural radiomic
information (e.g. energy, entropy, interquartile range) due to potential
overlap with functional radiomic information.

###  3.3 Methodology for extracting functional radiomics features via quantlet
basis approach

The starting point for the proposed scalar-on-functional quantile
classification approach is to compute the empirical quantiles. The quantile
distributions naturally account for heterogeneity in the image (arising from
lesional, peri-lesional, and healthy regions in the organ) and is able to
tackle different image sizes across samples. Figure 2 visualizes differences
in the quantile distributions between different lesion grades as well as
between cases vs controls. It is clear that the quantile distributions for the
lesion pixels are often flatter compared to the distributions for non-lesional
pixels. In particular, the non-lesional distribution has near-zero values and
higher values corresponding to left and right tails respectively, which
represents greater variability in the pixel distribution.

Figure 2:  Visualizations of lesion EPM images with segmentations for lesion
grades 2-5 and healthy control. Liver masks are shown in yellow; lesion masks
in red and surrounding peri-lesional areas are shaded in blue. Mean empirical
quantiles are shown along with 95% confidence bands for each lesion grade and
controls. For cases, quantiles are shown across all lesional (red) and peri-
lesional (blue) areas, and for controls, quantiles are shown across all liver
masks (yellow). Distributions of log-lesion volumes are also shown for each
lesion grade. Observed lesion volumes ranged from 0.17-0.950 cm3, with
approximate diameters of 0.70-12.2 cm.

First, we demonstrate how to derive the empirical quantiles from the EPM
images. Denote the $i$th image as $X_{i}$ comprising $m_{i}$ pixels
($i=1,\ldots,N$). We note that $X_{i}$ may include the full liver image for
subject level classification and a smaller image involving lesional and per-
lesional areas for lesion level classification. Further, denote the $m_{i}$
empirical quantiles for the $i$th sample as
$\\{Q_{i}(p_{1}),\ldots,Q_{i}(p_{m_{i}})\\}$, which are derived from $X_{i}$,
$i=1,\ldots,N,$ for $0<p_{j}<1$. Mathematically, the quantile function is
defined as $Q_{X}(p)=F^{-1}_{X}(p)=inf\\{x:F_{X}(x)\geq p\\},\mbox{ }0\leq
p\leq 1,$ where $inf(\cdot)$ denotes the infimum, and $F_{X}$ denotes the
cumulative distribution function, i.e. $p=F_{X}(x)=P(X\leq x)\in[0,1]$.
Further, denote the order statistics from the $i$th image as
$\\{X^{(1)}_{i},\ldots,X^{(m_{i})}_{i}\\}$ that are arranged in an increasing
order. We note that for $p\in[1/(m_{i}+1),m_{i}/(m_{i}+1)]$, the empirical
quantile function for the $i$th sample is computed as

$\displaystyle\hat{Q}_{i}(p)=(1-w_{i})X_{i}^{((m_{i}+1)p)}+w_{i}X_{i}^{((m_{i}+1)p+1)},$
(1)

where $w_{i}=X_{i}^{((m_{i}+1)p+1)}-X_{i}^{(m_{i}+1)p}$. The above expression
is obtained via interpolation, where $[x]$ is an integer $\leq x$. In
practice, we will focus our modeling on quantiles lying within
$[\delta,1-\delta]$ where $\delta=\max_{i\leq n}\\{1/(m_{i}+1)\\}$.

We note that the empirical quantile function can be computed on a grid
$\\{p_{1},\ldots,p_{K}\\}$ lying within $[\delta,1-\delta]$ using (1), where
the grid is standardized across images. Subsequently, the corresponding
empirical quantiles $\\{Q_{X}(p_{1}),\ldots,Q_{X}(p_{K})\\}$ could be
potentially used as covariate features for prediction. However, the quality of
the empirical quantiles computed over the standarized grid will depend on the
size of the images. While it may be possible to estimate these quantiles with
reasonable accuracy corresponding to larger images with several thousands of
pixels, the same may not be true for smaller images with fewer pixels or small
lesions. This is due to the fact that the order statistics
$\\{X^{(1)}_{i},\ldots,X^{(m_{i})}_{i}\\}$ are potentially spaced farther
apart for images with fewer pixels, that results in an inferior interpolation
in (1). Therefore, directly using the empirical quantile distribution as
functional predictors may not be desirable, particularly for lesion-level
analysis.

While useful, it may not be appealing to directly use the empirical quantiles
in the modeling approach due to the fact that their accuracy depends on the
image size and other factors such as lack of smoothness across neighboring
quantiles. In particular, the empirical quantile estimates may be sensitive to
the number of pixels in the image due to inadequate interpolation. A more
appealing approach is to propose a statistical smoothing method over the
empirical quantiles with the goal being to provide more accurate and robust
quantile-based features. As shown in [35], it is possible to find an empirical
set of basis functions $\\{\psi_{1},\ldots,\psi_{K}\\}$ on the quantile space
to approximate the empirical quantiles with extremely high accuracy (a
property known as near lossless). That is,

$\displaystyle
Q_{X_{i}}(p_{j})\approx\sum_{k=1}^{K}Q^{*}_{ik}\psi_{k}(p_{j}),$ (2)

where $Q^{*}_{ik}$ represent the basis coefficients for the $k$th basis
function for the $i$th sample, and $\\{\psi_{1}(p),\ldots,\psi_{K}(p)\\}$
represent a set of orthogonal quantlet basis functions. Under suitable choice
of basis functions, even a small to moderate number of basis functions in (2)
is sufficient to guarantee an accurate approximation, as shown in [35]. In
contrast to the empirical quantiles that is derived directly from the imaging
data, the functional quantlet approach uses a wavelet-based smoothing approach
that leads to information sharing across neighboring quantile levels.
Therefore, empirical quantiles are more susceptible to noise in the image and
the imaging size compared to the functional quantlet-based approach.

The basis coefficients in (2) are obtained in a straightforward manner as
follows. Let us denote $\bm{Q_{i}}=(Q_{i}(1),\ldots,Q_{i}(m_{i}))$, and
$\bm{\Psi}$ be a $K\times m_{i}$ with the $\bm{\Psi}(k,j)=\psi_{k}(p_{j})$.
Then the basis coefficients are given as
$\bm{Q}^{*}_{i}=(Q^{*}_{i1},\ldots,Q^{*}_{iK})=\bm{Q}_{i}\bm{\Psi}^{T}_{i}(\bm{\Psi}_{i}\bm{\Psi}^{T}_{i})^{-1}$,
where $K<<\min\\{m_{1},\ldots,m_{N}\\}$ in our settings involving a large
number of pixels in the image. We note $\bm{Q}^{*}_{i}$ contains virtually all
of the information in the raw data $\bm{Q}_{i}$ via the near lossless
property, but it is simultaneously more robust to noise and overfitting in the
empirical quantile distribution. Therefore, we use the smoothed quantile
function
$\hat{Q}_{X_{i}}(p_{j})=\sum_{k=1}^{K}\hat{Q}^{*}_{ik}\psi_{k}(p_{j})$ as the
functional predictor in our model in lieu of the empirical quantile
distribution, where $\hat{Q}^{*}$ denotes the estimated coefficients under
(2). More details on such coefficient estimation as well as the choice of $K$
are provided in the following subsection.

###  3.4 Choice of Quantlet basis

We adapt the approach in [35] for choosing the basis functions, and design
quantlets as a suitable set of basis functions that were shown to work well
for representing the quantile distribution. We create empirical quantlet bases
from a given dataset, serving as a robust representation of each subject’s
quantile function. The basis function is able to reconstruct the quantile
distribution with greater than $99.9\%$ accuracy using a moderate number
coefficient terms. This type of low dimensional feature extraction results in
massive dimension reduction for large images and yields a more faithful
reconstruction of the quantile distribution compared to alternate
representations such as Legendre-polynomials [11] that can result in poor tail
approximation for heavy tailed distributions [13]. These bases are employed in
our subsequent quantile functional regression model.

Consider the set of quantile functions $\mathcal{Q}$ such that
$\mathcal{Q}:=\\{Q:\int\big{(}Q^{2}(p)d\Pi(p)\big{)}^{1/2}<\infty\\}$, where
$\Pi(\cdot)$ is a uniform density with respect to the Lebesgue measure. Define
the first two basis functions to be a constant basis $\zeta_{1}(p)=1$ for
$p\in[0,1]$ and standard normal quantile function $\zeta_{2}(p)=\Phi^{-1}(p)$.
These orthogonal bases cover the entirety of Gaussian quantile functions,
where the initial coefficient signifies the mean and the second coefficient
represents the distribution’s variance. An overcomplete dictionary is created,
comprising these bases as well as numerous elements generated from Beta
cumulative density functions (CDFs). The additional basis functions are
constructed as
$\zeta_{j}(p)=P_{N^{\perp}}(\frac{F_{\theta_{k}}(p)-\mu_{\theta_{k}}}{\sigma_{\theta_{k}}}),j\geq
3,$ where $F_{\theta_{k}}$ is the cumulative distribution function (CDF) of a
Beta$(\theta_{k})$ distribution for some positive parameters
$\theta_{k}=\\{a_{k},b_{k}\\}$,
$\mu_{\theta_{k}}=\int_{0}^{1}F_{\theta_{k}}(u)du$, and
$\sigma^{2}_{\theta_{k}}=\int_{0}^{1}(F_{\theta_{k}}(u)-\mu_{\theta_{k}})^{2}du$
are the centered and scaled values of these distributions for standardization,
respectively, and $P_{N^{\perp}}$ indicates the projection operator onto the
orthogonal complement to the Gaussian basis elements $\zeta_{1}$ and
$\zeta_{2}$. In other words,
$P_{N^{\perp}}(f(p))=f(p)-\zeta_{1}(p)\int_{0}^{1}f(u)\zeta_{1}(u)du-\zeta_{2}(p)\int_{0}^{1}f(u)\zeta_{2}(u)du$.

Put together, the set
$\mathcal{D}^{O}=\\{\zeta_{1}\cup\zeta_{2}\\}\cup\\{\zeta_{k}:\theta_{k}\in\mathcal{\theta}\\}_{k=1}^{K_{0}}$
represents an overcomplete dictionary, where $\mathcal{\theta}$ represents the
parameter space of interest. In practice, to fix the number of dictionary
elements, we choose a grid on the parameter space to obtain
$\mathcal{\theta}=\\{\theta_{k}=(a_{k},b_{k})\\}_{k=3}^{K_{0}}$ by uniformly
sampling on $\mathcal{\theta}\subset(0,J)^{2}$ for some sufficiently large
$J$, and and choosing $K_{0}$ to be a large integer. More details on how to
choose the parameter space is provided in [35], along with additional
theoretical results that justify the choice of a large dictionary of Beta CDF
to represent any quantile function whose first derivative is absolutely
continuous.

Once this overcomplete basis set has been constructed, a sparse subset of
basis functions are chosen that can represent the quantile functions across
all samples. This is made possible by first fitting a Lasso regression [31] to
select a sparse set of wavelets for each subject separately, and subsequently
taking the union across these sparse sets of basis elements across all
subjects, to obtain a unified set of dictionary elements $\mathcal{D}$.
Subsequently, a leave-one-out concordance correlation coefficient (LOOCCC) is
proposed that allows one to reduce the number of basis elements in
$\mathcal{D}$ even further, with minimal loss of information. This property is
known as near lossless property. This goodness of fit criteria for the $i$th
subject is given as
$\rho_{(i)}=\frac{Cov(Q_{i}(\cdot),\sum_{k\in\mathcal{D}_{(i)}^{U}}\zeta_{k}(\cdot)Q^{U}_{ik})}{Var(Q_{i}(\cdot))+Var(\sum_{k\in\mathcal{D}_{(i)}^{U}}\zeta_{k}(\cdot)Q^{U}_{ik})}$,
where $\mathcal{D}_{(i)}^{U}$ is the set of the sparse basis elements obtained
under the Lasso regression, combined across all but the $i$-th subject.
Subsequently one can use $\rho^{0}=\min\\{\rho_{(i)}\\}_{i}$ as an overall
goodness of fit measure that can be used to select a small to moderate number
of basis elements for our downstream prediction model, with no more than
minimal loss of information For example, in the Liver EPM dataset, we find
that $13$ basis elements are able to preserve a concordance of at least 0.999,
when reconstructing the quantile distributions from the entire liver mask.
Once a reduced set of basis elements are obtained, subsequent Gram-Schimdt
orthogonalization steps are used followed by a wavelet based denoising step
[6].

###  3.5 Classification based on extracted radiomics features

####  3.5.1 Proposed Model

We generalize the classical logistic regression models [25] to include
functional predictors in the form of smoothed quantile distributions. We use
the derived functional quantlet features, along with the extracted structural
radiomics features and demographic features for risk prediction and early
diagnosis. Let $y_{i}\sim Ber(p_{i}(Q_{x}))$ denote the binary outcome of
interest, which follows a Bernoulli distribution with probability
$p_{i}(Q_{x})$ that depends on the quantile distribution of the image.
Further, let ${\bf y}=(y_{1},\ldots,y_{N})$ denote the vector of binary
outcomes, ${\bf w}$ denotes the demographic features, and ${\bf z}$ denotes
the structural radiomics features derived from the arterial phase of the
T1w-MRI scans.  Let $\mathcal{X}_{i}$ correspond to the full liver image for
Aim 1A, and let it correspond to pre-defined ROIs encompassing the lesion for
the cases in Aim 1B and for the lesion level analyses in Aims 2-3. When lesion
segmentation is available, denote the lesion mask for the $i$th image as
$\mathcal{L}_{i}$, so that $\mathcal{L}_{i}^{c}=X_{i}\setminus\mathcal{L}_{i}$
denotes the peri-lesional areas adjacent to the lesion. The likelihood is
$L({\bf
y}\mid\mathcal{X})=\prod_{i=1}^{N}p^{y_{i}}_{i}(Q_{x})(1-p_{i}(Q_{x}))^{1-y_{i}}$,
where

$\displaystyle\frac{p_{i}(Q_{x})}{1-p_{i}(Q_{x})}=\alpha+({\bf w}_{i}^{T},{\bf
z}_{i}^{T})\bm{\zeta}+\int_{0}^{1}\hat{Q}_{\mathcal{X}_{i}}(p)\beta_{\mathcal{X}}(p)dp$
$\displaystyle\approx\alpha+({\bf w}_{i}^{T},{\bf
z}_{i}^{T})\bm{\zeta}+\sum_{k=1}^{K}\hat{Q}^{*}_{\mathcal{L}_{i}k}\gamma_{k}+\sum_{k=1}^{K}\hat{Q}^{*}_{\mathcal{L}_{i}^{c},k}\eta_{k},$
(3)

where
$\int_{0}^{1}\hat{Q}_{\mathcal{X}_{i}}(p)\beta_{\mathcal{X}}(p)dp=\int\hat{Q}_{\mathcal{L}_{i}}(p)\beta_{\mathcal{L}}(p)dp+\int\hat{Q}_{\mathcal{L}^{c}_{i}}(p)\beta_{\mathcal{L}^{c}}(p)dp\approx\sum_{k=1}^{K}\hat{Q}^{*}_{\mathcal{L}_{i}k}\gamma_{k}+\sum_{k=1}^{K}\hat{Q}^{*}_{\mathcal{L}_{i}^{c},k}\eta_{k}$,
$\beta(p)$ denotes the functional regression coefficient that can be
represented as the basis expansion
$\beta(p)\approx\sum_{k=1}^{K}\eta_{k}\psi_{k}(p)$, $Q^{*}$ are defined in
(2), ${\bf w}$ denotes the demographic features with corresponding effects
denoted by ${\bm{\zeta}}$, and ${\bf z}$ denotes the structural radiomics
features derived from the arterial phase of the T1w-MRI scans with
corresponding effects $\bm{\gamma}$. Model (3) accounts for heterogeneity by
assigning separate effects for the lesional and peri-lesional areas, and via
smoothed quantile curves capturing rich distributional information across all
pixels. We the following penalized likelihood criteria to estimate the model
parameters in (3):

$\displaystyle\hat{\bm{\theta}}=\mbox{arg min}_{\bm{\theta}}-\log(L({\bf
y}\mid\mathcal{X},\bm{\theta})+\lambda|\bm{\theta}|_{1},$ (4)

where $\bm{\theta}=(\alpha,\bm{\zeta,\gamma,\eta})$ denotes the collection of
model parameters, $\lambda$ denotes the penalty parameter for the $L_{1}$
penalty, with a larger value indicating greater sparsity and vice-versa. The
imposed $L_{1}$ penalty down-weights the unimportant coefficients to zero,
thus reducing overfitting and resulting in model parsimony. The penalized
logistic regression model in (3)-(4) was implemented in R using the glmnet
package. The regularization parameter $\lambda$ was selected by using 10-fold
cross-validation on the entire dataset, and this choice of $\lambda$ is
subsequently used in all model fitting to report the leave-out-one cross-
validation prediction errors. We denote the proposed model in (3)-(4) as
scalar-on-functional quantile classification model that is used to address the
following analysis goals.

####  3.5.2 Description of Analysis Goals

We focus on three distinct aims with strong translational implications. The
first goal (Aim 1: Diagnostic Aim) involves differentiating between cases and
controls using the Liver EPM images at both pre-diagnostic and diagnostic
visits using a cross-sectional analysis. We note that all the control subjects
had one visit, whereas individuals with lesions had two visits. This resulted
in a classification analysis with 97 control subjects and 91 images with
lesions. Aim 1A uses the whole liver EPM images without knowledge of the
lesion masks when computing the functional quantlet features to predict cases
versus controls. This set-up is challenging since the information from the
lesional tissue is somewhat masked by the information contained in the healthy
tissues that encompasses most of the liver, especially given the relatively
small lesion sizes. The resulting analysis contains limited information about
the lesion characteristics that potentially yields reduced power to
differentiate cases versus controls. We further extend the analysis to include
information about the lesion masks for cases in Aim 1B. Aim 1B is expected to
produce superior classification accuracy since the quantile distributions from
the cases arise from lesional and peri-lesional areas only, which provides a
greater differentiating power between cases and controls.  The second goal
(Aim 2: Risk Stratification) focuses on lesion level classification. In
particular, we are interested in differentiating between aggressive tumors vs
mild tumors in a cross-sectional manner by combining information across both
visits but without explicitly accounting for longitudinal dependence within
subjects. A tumor is labeled as aggressive if it has grade higher than 3, and
mild if the lesion grade is 3 or lower. This results in a sample size of 84
for this analysis, with 39 tumors having grades higher than 3, and the rest
having lower grades.

Transitioning beyond Aims 1-2 that are designed primarily for demonstrating
the practicality of the proposed functional quantlet features, we now delve
into our ultimate and pivotal goal of early detection. Aim 3 specifically
concentrates on forecasting tumor progression during the diagnostic visit
using imaging data solely from the pre-diagnostic visit. The primary goal is
to discern which lesions are likely to remain stable over time compared to
those anticipated to deteriorate, indicated by an escalation in tumor grade
between visits. This analysis is conducted at the lesion level, involving 35
samples. The modest sample size presents a challenge for most models to
effectively predict lesion progression, particularly when considering the
large size of the images.

In order to include information from the adjoining areas surrounding the
lesions that can potentially be impacted by the lesion in the future, we
construct a sphere-shaped region of interest (ROI) around the lesion that
encompasses the entire lesion and is approximately 2 times the volume of the
corresponding lesion. These areas, referred to as peri-lesional areas, are
likely to be significantly influenced by tumor progression but remain
unaffected for stable tumors that exhibit minimal changes over time.
Consequently, they are anticipated to contain valuable information crucial to
our objective of early detection and are thus integral to our analysis.
Notably, the inclusion of these adjacent peri-lesional areas resulted in
considerable improvements in classification accuracy for Aim 3 compared to a
more restrictive analysis focused strictly on the lesion mask (results not
presented here).

####  3.5.3 Benchmark approaches for comparison

A leave-out-one cross-validation (LOOCV) strategy was used to evaluate the
performance of the proposed approach versus various competing models that
differed in terms of the sets of predictors used for classification under a
penalized logistic regression approach. In particular, we compared with
classification analyses involving: (i) demographic covariates only; (ii)
summary measures such as mean or median of the pixel-level EPM values along
with demographics; (iii) empirical quantiles along with demographics; (iv)
summary structural radiomics features with demographics; and (v) the features
considered in (ii)-(iii) along with structural radiomics features. We consider
various numbers of basis coefficients corresponding to classification analysis
based on quantiles.

Several metrics were used to evaluate performance. These measures included:
(a) sensitivity that is defined as the power to detect true cases in Aim 1,
the power to detect tumors with higher grades in Aim 2, and the power to
detect aggressive tumors that longitudinally progresses to higher grades in
Aim 3; (b) the corresponding specificity rates that is equivalent to 1- false
discovery rate; (c) the F1 score that is the harmonic mean of sensitivity and
specificity; and (d) accuracy that is defined as the proportion of correctly
identified cases and controls. In our overall analysis, we fit the models with
the extracted radiomics features as well as using principal components derived
from these features that are designed to remove collinearity. However, the
latter PCA based analysis were not included since they did not produce
improvements.

##  4 Results

The results for Aims 1-3 are presented in Tables 1-3, respectively. We report
the results over varying number of basis coefficients included in the
classification model, where these coefficients are included based on their
ranking (in a decreasing order of importance) in terms of their LOOCCC score.
From Table 1, it is clear that the best classification accuracy is obtained
when using the functional quantlet features along with the structural
radiomics features and demographics. For the case vs classification model, we
observe that combining the functional quantlet features with the structural
radiomics features results in considerable gains in sensitivity along with
comparable or improved specificity compared to using either type of features
alone. Moreover, classification based on quantlet features along with
structural radiomics features results in improvements in sensitivity compared
to classification based on empirical quantiles as well as summary mean(median)
EPM values, when they are used in isolation or coupled with demographics.

| Cases vs Control Diagnosis
---|---
| Aim 1A | Aim 1B
Features | Sens | Spec | F1 | Acc | Sens | Spec | F1 | Acc
D | 0.07 | 0.94 | 0.12 | 0.57 | 0.07 | 0.94 | 0.12 | 0.57
Mean | 0.36 | 0.76 | 0.43 | 0.59 | 0.51 | 0.90 | 0.62 | 0.73
Median | 0.29 | 0.80 | 0.38 | 0.58 | 0.63 | 0.93 | 0.73 | 0.80
EQ(10) | 0.26 | 0.95 | 0.40 | 0.65 | 0.63 | 0.98 | 0.76 | 0.83
EQ(20) | 0.26 | 0.93 | 0.39 | 0.64 | 0.64 | 0.99 | 0.77 | 0.84
EQ(30) | 0.29 | 0.93 | 0.42 | 0.65 | 0.68 | 1.00 | 0.81 | 0.86
EQ(50) | 0.29 | 0.92 | 0.42 | 0.65 | 0.70 | 1.00 | 0.83 | 0.88
Q(10) | 0.38 | 0.79 | 0.45 | 0.61 | 0.89 | 0.93 | 0.90 | 0.91
Q(20) | 0.42 | 0.76 | 0.48 | 0.61 | 0.89 | 0.92 | 0.89 | 0.90
Q(30) | 0.44 | 0.74 | 0.50 | 0.61 | 0.89 | 0.93 | 0.90 | 0.91
R | 0.75 | 0.94 | 0.82 | 0.86 | 0.75 | 0.94 | 0.82 | 0.86
Mean+R | 0.75 | 0.94 | 0.82 | 0.86 | 0.82 | 0.96 | 0.87 | 0.90
Median+R | 0.75 | 0.94 | 0.82 | 0.86 | 0.82 | 0.96 | 0.87 | 0.90
EQ(10)+R | 0.71 | 0.97 | 0.81 | 0.86 | 0.76 | 0.98 | 0.85 | 0.89
EQ(20)+R | 0.71 | 0.97 | 0.81 | 0.86 | 0.76 | 0.98 | 0.85 | 0.89
EQ(30)+R | 0.71 | 0.97 | 0.81 | 0.86 | 0.75 | 0.99 | 0.85 | 0.89
EQ(50)+R | 0.71 | 0.97 | 0.81 | 0.86 | 0.79 | 0.99 | 0.88 | 0.90
Q(10)+R | 0.78 | 0.95 | 0.84 | 0.88 | 0.88 | 0.98 | 0.92 | 0.93
Q(20)+R | 0.74 | 0.97 | 0.83 | 0.87 | 0.88 | 0.98 | 0.92 | 0.93
Q(30)+R | 0.83 | 0.93 | 0.86 | 0.89 | 0.86 | 0.98 | 0.91 | 0.93
Table 1: Aim 1 Results: Leave-out-one CV results for case/control classification. In Aim 1A, EPM features are taken from liver masks for control and cases, whereas for Aim 1B, EPM features for the cases are taken from their lesion masks. Different predictors were considered including Mean and Median EPM, EPM empirical quantiles (EQ), EPM quantlets (Q), and PyRadiomics features (R). We considered varying number of quantiles and 70 PyRadiomics features (R). | Risk Stratification
---|---
Features | Sens | Spec | F1 | Acc
D | 0.55 | 0.58 | 0.58 | 0.56
Mean | 0.63 | 0.27 | 0.56 | 0.47
Median | 0.65 | 0.30 | 0.58 | 0.49
EQ(10) | 0.63 | 0.36 | 0.58 | 0.51
EQ(30) | 0.65 | 0.27 | 0.58 | 0.48
EQ(50) | 0.60 | 0.30 | 0.55 | 0.47
Q(10) | 0.68 | 0.36 | 0.61 | 0.53
Q(30) | 0.60 | 0.52 | 0.60 | 0.56
Q(50) | 0.83 | 0.58 | 0.76 | 0.71
R | 0.80 | 0.64 | 0.76 | 0.73
Mean+R | 0.80 | 0.64 | 0.76 | 0.73
Median+R | 0.80 | 0.64 | 0.76 | 0.73
EQ(10)+R | 0.78 | 0.67 | 0.76 | 0.73
EQ(30)+R | 0.78 | 0.67 | 0.76 | 0.73
EQ(50)+R | 0.78 | 0.67 | 0.76 | 0.73
Q(10)+R | 0.83 | 0.64 | 0.77 | 0.74
Q(30)+R | 0.85 | 0.64 | 0.79 | 0.75
Q(50)+R | 0.80 | 0.76 | 0.80 | 0.78
Table 2: Aim 2 Results: Leave-out-one CV results for cross-sectional
classification of lesion grade.

Table 1 shows that while the sensitivity is greater than 80% and specificity
is higher than 90% under both Aims 1A-1B, the performance is superior for Aim
1B. The improved performance for Aim 1B is due to the fact that the quantile
distributions arising from the lesional and peri-lesional areas were
considerably more different than the quantile distributions arising from
healthy liver parenchyma. In contrast, the entire liver image was used for Aim
1A without any lesion level information that may have diluted the
discriminative ability of the methods. For Aim 2 involving lesion grade
classification (see Table 2), similar trends hold. In particular, combining
the functional quantlet features with the structural radiomics features
results in sensitivity of 80% that is higher than using each type of features
alone in the model, and specificity of 76% with 50 quantlets. In contrast,
using the empirical quantiles, as well as summary measures (mean,median),
results in lower sensitivity and specificity, when used alone or in
conjunction with structural radiomics features. This points to the limited
utility of these features in risk stratification.

| Progression (lesion) | Progression (ROI)
---|---|---
Features | Sens | Spec | F1 | Acc | Sens | Spec | F1 | Acc
D | 0.18 | 0.90 | 0.27 | 0.65 | 0.18 | 0.90 | 0.27 | 0.65
Mean | 0 | 0.83 | 0 | 0.54 | 0.18 | 0.85 | 0.25 | 0.61
Median | 0.10 | 0.89 | 0.15 | 0.61 | 0.27 | 0.80 | 0.33 | 0.61
EQ(10) | 0 | 0.94 | 0 | 0.61 | 0.46 | 0.70 | 0.46 | 0.61
EQ(20) | 0.10 | 0.94 | 0.17 | 0.64 | 0.46 | 0.80 | 0.50 | 0.68
EQ(50) | 0 | 0.89 | 0 | 0.57 | 0.46 | 0.80 | 0.50 | 0.68
Q(10) | 0 | 0.89 | 0 | 0.57 | 0.46 | 1.00 | 0.63 | 0.81
Q(20) | 0.20 | 0.83 | 0.27 | 0.61 | 0.46 | 1.00 | 0.63 | 0.81
Q(50) | 0.50 | 0.94 | 0.62 | 0.79 | 0.27 | 0.95 | 0.40 | 0.71
R | 0.60 | 0.94 | 0.71 | 0.82 | 0.64 | 0.90 | 0.70 | 0.81
Mean+R | 0.60 | 0.94 | 0.71 | 0.82 | 0.64 | 0.80 | 0.64 | 0.74
Median+R | 0.60 | 0.94 | 0.71 | 0.82 | 0.64 | 0.90 | 0.70 | 0.81
EQ(10)+R | 0.60 | 1 | 0.75 | 0.86 | 0.64 | 0.90 | 0.70 | 0.81
EQ(20)+R | 0.60 | 1 | 0.75 | 0.86 | 0.64 | 0.90 | 0.70 | 0.81
Q(10)+R | 0.80 | 0.94 | 0.84 | 0.90 | 0.82 | 0.85 | 0.78 | 0.84
Q(20)+R | 0.80 | 0.94 | 0.84 | 0.90 | 0.82 | 0.85 | 0.78 | 0.84
Table 3: Aim 3 Results: Leave-out-one CV results for change in lesion grade
across visits ($>1$ vs. $\leq 1$). Two separate analysis are performed, one
exclusively involving the lesion voxels, and the other involving the ROI that
contains additional peri-lesional tissues surrounding the lesion.

For predicting tumor progression (Aim 3), we report results from two separate
analysis in Table 3. In the first analysis the prediction is based on pixels
lying within the lesion exclusively (denoted as lesional analysis), whereas
for the second analysis, the prediction is based on pre-defined ROIs that
contain some peri-lesional areas surrounding the lesion (denoted as ROI
analysis). We find that the best results are obtained under the proposed
approach combining the functional quantlet and structural radiomics features
derived from the lesions directly, without including peri-lesional areas. In
particular, the analysis focusing on the lesional areas exclusively produces
comparable sensitivity (80%) but much higher specificity (94%) under the
proposed approach when compared to the corresponding ROI level analysis
involving additional peri-lesional areas. We therefore conclude that including
peri-lesional areas surrounding the lesion leads to loss of specificity (or
equivalently an increase in the false discovery rate) without providing any
obvious benefits in terms of an increase in sensitivity. Moreover, alternate
analysis involving empirical quantiles and summary measures produce
considerably lesser sensitivity and comparable specificity under both types of
analysis (lesional analysis as well as ROI analysis). Finally, the F-1 score
is markedly higher for the lesional analysis under the proposed approach.
Coupled with high sensitivity levels, these results highlight the prognostic
value of the imaging biomarkers.

Overall, for all the three aims, a combination of the functional quantlet
features and structural radiomics features produced the best results with
increased sensitivity and comparable specificity that translated to improved
accuracy and F-1 scores. The improvements under the proposed smoothed quantile
approach over empirical quantiles is more noticeable for Aims 2 and 3 where
the lesion mask information was incorporated. Moreover, combining the
functional quantlet features with structural radiomics produces greater
improvements in classification compared to the combination of empirical
quantiles and structural radiomics features. This suggests that the smoothed
quantile distribution can leverage informative patterns from the lesion for
classification when the lesion location is specified, and that they are more
amenable for integrative analysis in combination with structural imaging
features. For Aim 3, the sharp increase in sensitivity for the lesional
analysis points to the ability of the proposed approach (combining functional
quantlet and structural radiomic features) to detect lesions that are likely
to grow aggressively over time. In contrast, the poor sensitivity of the
routinely used summary measures (mean, median) is not desirable and points to
the limits of existing analysis approaches in terms of accommodating
heterogeneity across the image.

##  5 Discussion

In this paper, we proposed a novel DDA based approach for image-based
prediction that accounts for image heterogeneity and incorporates information
from the pixel-level quantile distribution as functional covariates, along
with summary structural radiomics features. Our extensive numerical analysis
revealed that the functional quantlet approach yielded improvements in terms
of out-of-sample classification, owing to the ability to capture minute
distributional differences often ignored by classical approaches. The
overarching framework developed here for image based risk prediction and early
diagnosis is expected to have considerable impact in cancer imaging. Some
additional considerations are elaborated below.

Statistical Consideration and Potential Pitfalls: For the lesion-level
classification, we attempted to utilize radiomic information only within or
directly surrounding each lesion for cross-sectional and longitudinal
prediction. However, one challenge that arises is the vast amount of
heterogeneity in the size (i.e., number of pixels) comprising the different
lesions. In this sense, our analysis was limited with respect to smaller
lesions which did not contain as much information. We note potential bias
arising from the fact that lower severity lesions tended to have smaller
sizes. Furthermore, these smaller lesions limited the scope of our feature
space both for the empirical quantiles, and the quantlet coefficients, which
were estimated over a unified grid, to limit cost and avoid additional issues.

In addition to the lesion size differences, we note that the choice of ROIs
representing the peri-lesional areas may impact the final results. In our
analysis, we selected ROIs as spherical regions centered around the lesion
centroid, with approximately 2 times the volume of the lesion. However, the
sensitivity of the approach to the choice of size and shape of these ROIs need
to be investigated further. Moreover, the analysis in Aim 3 demonstrated that
the inclusion of surrounding peri-lesional tissues is not always guaranteed to
produce improvements compared to an analysis focused exclusively on lesional
pixels. Ultimately, the performance may depend on additional considerations
such as sample size. Finally, the recommended number of quantlets to be
included in the model should be such that the proposed method achieves at
least 80% sensitivity and specificity in a validation sample.

Computational Scalability: The computational scalability of functional
quantlet classification relies largely on the process of estimating quantlet
coefficients. The computational speed of the quantlet extraction process
depends on factors such as the variability in image sizes. When the image
sizes vary widely, the quality of the empirical quantiles computed from the
smaller images may be compromised resulting in computation challenges related
to issues with non-uniform sampling, Gram-Schmidt orthogonalization, and
potential memory issues in the basis dictionary construction. To get around
these challenges, we computed quantlets based on $512$ empirical quantiles for
the case-control classification (entire liver) and $128$ empirical quantiles
for the lesion-level classifications. While this downsampling may result is
minimal information loss, it facilitates practical implementation and
stability.

Future Directions: The framework developed in this article can be generalized
in several different directions. For example, one could evaluate whether
longitudinal changes in the EPM images are associated with tumor progression.
The proposed framework can also be used to predict response to therapy in
future visits, in clinical trials. Future analysis could also be enriched by
including additional non-imaging biomarkers such as blood based biomarkers
frequently used in HCC. Finally, the proposed DDA methodology can be broadly
applied to other types of imaging modalities and disease areas.

## References

  * [1] Hugo JWL Aerts, Emmanuel Rios Velazquez, Ralph TH Leijenaar, Chintan Parmar, Patrick Grossmann, Sara Carvalho, Johan Bussink, René Monshouwer, Benjamin Haibe-Kains and Derek Rietveld  “Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach”  In _Nat Comm_ 5 Nature Publishing Group UK London, 2014, pp. 4006
  * [2] Nicole H Augustin, Calum Mattocks, Julian J Faraway, Sonja Greven and Andy R Ness  “Modelling a response as a function of high-frequency count data: the association between physical activity and fat mass”  In _Stat Meth in Med Res_ 26.5 Sage Publications Sage UK: London, England, 2017, pp. 2210–2226
  * [3] Madhuri Avula, Narasimha Prasad Lakkakula and Murali Prasad Raja  “Bone Cancer Detection from MRI Scan Imagery Using Mean Pixel Intensity”  In _2014 8th Asia Modelling Symposium_ , 2014, pp. 141–146  DOI: 10.1109/AMS.2014.36
  * [4] Victoria Chernyak, Kathryn J Fowler, Aya Kamaya, Ania Z Kielar, Khaled M Elsayes, Mustafa R Bashir, Yuko Kono, Richard K Do, Donald G Mitchell and Amit G Singal  “Liver Imaging Reporting and Data System (LI-RADS) version 2018: imaging of hepatocellular carcinoma in at-risk patients”  In _Radiology_ 289.3 Radiological Society of North America, 2018, pp. 816–830
  * [5] Ibiayi Dagogo-Jack and Alice T Shaw  “Tumour heterogeneity and resistance to cancer therapies”  In _Nat Rev Clin Onc_ 15.2 Nature Publishing Group UK London, 2018, pp. 81–94
  * [6] David L Donoho, Iain M Johnstone, Gérard Kerkyacharian and Dominique Picard  “Wavelet shrinkage: asymptopia?”  In _J of the Royal Stat Soc: Series B_ 57.2 Wiley Online Library, 1995, pp. 301–337
  * [7] Dorothea Dumuid, Željko Pedišić, Javier Palarea-Albaladejo, Josep Antoni Martín-Fernández, Karel Hron and Timothy Olds  “Compositional data analysis in time-use epidemiology: what, why, how”  In _Int J of Env Res and Pub Health_ 17.7 MDPI, 2020, pp. 2220
  * [8] Xiangnan Feng, Tengfei Li, Xinyuan Song and Hongtu Zhu  “Bayesian scalar on image regression with nonignorable nonresponse”  In _J Am Stat Assoc_ 115.532 Taylor & Francis, 2020, pp. 1574–1597
  * [9] Karl J Friston, Chris D Frith, Robert Turner and Richard SJ Frackowiak  “Characterizing evoked hemodynamics with fMRI”  In _Neuroimage_ 2.2 Elsevier, 1995, pp. 157–165
  * [10] Laya Ghodrati and Victor M Panaretos  “Distribution-on-distribution regression via optimal transport maps”  In _Biometrika_ 109.4 Oxford University Press, 2022, pp. 957–974
  * [11] Rahul Ghosal, Vijay R Varma, Dmitri Volfson, Inbar Hillel, Jacek Urbanek, Jeffrey M Hausdorff, Amber Watts and Vadim Zipunnikov  “Distributional data analysis via quantile functions and its application to modeling digital biomarkers of gait in Alzheimer’s disease”  In _Biostatistics_ 24.3 Oxford University Press, 2023, pp. 539–561
  * [12] Rajarshi Guhaniyogi, Shaan Qamar and David B Dunson  “Bayesian tensor regression”  In _J Mach Learn Res_ 18.1 JMLR. org, 2017, pp. 2733–2763
  * [13] Jonathan RM Hosking  “L-moments: analysis and estimation of distributions using linear combinations of order statistics”  In _J of the Royal Stat Soc, Series B_ 52.1 Oxford University Press, 1990, pp. 105–124
  * [14] The MathWorks Inc.  “MATLAB version: 9.13.0 (R2022b)”  Natick, Massachusetts, United States: The MathWorks Inc., 2022  URL: https://www.mathworks.com
  * [15] Suprateek Kundu, Alec Reinhardt, Serena Song, Joo Han, M Lawson Meadows, Bruce Crosson and Venkatagiri Krishnamurthy  “Bayesian longitudinal tensor response regression for modeling neuroplasticity”  In _Human Brain Mapping_ Wiley Online Library, 2023
  * [16] Philippe Lambin, Ralph TH Leijenaar, Timo M Deist, Jurgen Peerlings, Evelyn EC De Jong, Janita Van Timmeren, Sebastian Sanduleanu, Ruben THM Larue, Aniek JG Even and Arthur Jochems  “Radiomics: the bridge between medical imaging and personalized medicine”  In _Nat Rev Clin Onc_ 14.12 Nature Publishing Group UK London, 2017, pp. 749–762
  * [17] Philippe Lambin, Emmanuel Rios-Velazquez, Ralph Leijenaar, Sara Carvalho, Ruud GPM Van Stiphout, Patrick Granton, Catharina ML Zegers, Robert Gillies, Ronald Boellard and André Dekker  “Radiomics: extracting more information from medical images using advanced feature analysis”  In _Euro. J. of Cancer_ 48.4 Elsevier, 2012, pp. 441–446
  * [18] Ruben THM Larue, Gilles Defraene, Dirk De Ruysscher, Philippe Lambin and Wouter Van Elmpt  “Quantitative radiomics studies for tissue characterization: a review of technology and methodological procedures”  In _The British J of Rad_ 90.1070 The British Institute of Radiology., 2017, pp. 20160665
  * [19] Kuo-Jung Lee, Galin L Jones, Brian S Caffo and Susan Spear Bassett  “Spatial Bayesian variable selection models on functional magnetic resonance imaging time-series data”  In _Bayesian Analysis_ 9.3 NIH Public Access, 2014, pp. 699
  * [20] Xin Ma, Suprateek Kundu and Alzheimer’s Disease Neuroimaging Initiative  “Multi-task learning with high-dimensional noisy images”  In _J Am Stat Assoc_ Taylor & Francis, 2022, pp. 1–14
  * [21] Nadim Mahmud, Maarouf A Hoteit and David S Goldberg  “Risk Factors and Center-Level Variation in Hepatocellular Carcinoma Under-Staging for Liver Transplantation”  In _Liver Trans_ 26.8 Wiley Online Library, 2020, pp. 977–988
  * [22] M Matabuena and A Petersen  “Distributional data analysis with accelerometer data in a nhanes database with nonparametric survey regression models.”  In _arXiv:2104.01165_ , 2021
  * [23] Jeffrey S Morris  “Functional regression”  In _Ann. Rev Stat and Its Appl_ 2 Annual Reviews, 2015, pp. 321–359
  * [24] Newsha Nikzad, D Fuentes, Millicent Roach, Tasadduk Chowdhury, Matthew Cagley, Mohamed Badawy, Manal M Hassan, K Elsayes, Laura Beretta and Eugene J Koay  “Characterization of the Imaging Signature of Hepatocellular Carcinoma with Enhancement Pattern Mapping”  In _HEPATOLOGY_ 76, 2022, pp. S1356–S1356  WILEY 111 RIVER ST, NJ USA
  * [25] Xianli Pan and Yitian Xu  “A Safe Feature Elimination Rule for $L\\_${$1$} L 1-Regularized Logistic Regression”  In _IEEE Trans on Pat Anal and Mach Intel_ 44.9 IEEE, 2021, pp. 4544–4554
  * [26] Peter C Park, Gye W Choi, Mohamed M., Dalia Elganainy, Danyal A Smani, John Tomich, Ray Samaniego, Jingfei Ma, Eric P Tamm and Sam Beddar  “Enhancement pattern mapping technique for improving contrast-to-noise ratios and detectability of hepatobiliary tumors on multiphase computed tomography”  In _Med Phys_ 47.1 Wiley Online Library, 2020, pp. 64–74
  * [27] James O Ramsay and CJ1125714 Dalzell  “Some tools for functional data analysis”  In _J. of the Royal Stat Soc Series B_ 53.3 Oxford University Press, 1991, pp. 539–561
  * [28] Philip T Reiss, Lan Huo, Yihong Zhao, Clare Kelly and R Todd Ogden  “Wavelet-domain regression and predictive inference in psychiatric neuroimaging”  In _The Ann of Appl Stat_ 9.2 NIH Public Access, 2015, pp. 1076
  * [29] Michael Smith and Ludwig Fahrmeir  “Spatial Bayesian variable selection with application to functional magnetic resonance imaging”  In _J. Am. Stat. Assoc._ 102.478 Taylor & Francis, 2007, pp. 417–431
  * [30] Renáta Talská, Karel Hron and Tomáš Matys Grygar  “Compositional scalar-on-function regression with application to sediment particle size distributions”  In _Math Geosciences_ 53.7 Springer, 2021, pp. 1667–1695
  * [31] Robert Tibshirani  “Regression shrinkage and selection via the lasso”  In _J of the Roy Stat Soc: Series B_ 58.1 Wiley Online Library, 1996, pp. 267–288
  * [32] Kristina Tzartzeva, Joseph Obi, Nicole E Rich, Neehar D Parikh, Jorge A Marrero, Adam Yopp, Akbar K Waljee and Amit G Singal  “Surveillance imaging and alpha fetoprotein for early detection of hepatocellular carcinoma in patients with cirrhosis: a meta-analysis”  In _Gastroenterology_ 154.6 Elsevier, 2018, pp. 1706–1718
  * [33] Joost JM Van Griethuysen, Andriy Fedorov, Chintan Parmar, Ahmed Hosny, Nicole Aucoin, Vivek Narayan, Regina GH Beets-Tan, Jean-Christophe Fillion-Robin, Steve Pieper and Hugo JWL Aerts  “Computational radiomics system to decode the radiographic phenotype”  In _Can Res_ 77.21 AACR, 2017, pp. e104–e107
  * [34] Xiao Wang, Hongtu Zhu and Alzheimer’s Disease Neuroimaging Initiative  “Generalized scalar-on-image regression models via total variation”  In _J Am Stat Assoc_ 112.519 Taylor & Francis, 2017, pp. 1156–1168
  * [35] Hojin Yang, Veerabhadran Baladandayuthapani, Arvind UK Rao and Jeffrey S Morris  “Quantile function on scalar regression analysis for distributional data”  In _J Am Stat Assoc_ 115.529 Taylor & Francis, 2020, pp. 90–106
  * [36] Hasan Çağrı Yıldırım, Gozde Kavgaci, Elvin Chalabiyev and Omer Dizdar  “Advances in the Early Detection of Hepatobiliary Cancers”  In _Cancers_ 15.15 MDPI, 2023, pp. 3880
  * [37] Vadim Zipunnikov, Brian Caffo, David M Yousem, Christos Davatzikos, Brian S Schwartz and Ciprian Crainiceanu  “Functional principal component model for high-dimensional brain imaging”  In _NeuroImage_ 58 Elsevier, 2011, pp. 772–784