#  The Problem of Marginality in Model Reductions of Turbulence

C. Gillot  CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France  G. Dif-Pradalier
<EMAIL_ADDRESS>CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France
Y. Sarazin  CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France  C. Bourdelle
CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France  Y. Camenen  Aix-Marseille
Université, CNRS PIIM, UMR 7345 Marseille, France  J. Citrin  DIFFER—Dutch
Institute for Fundamental Energy Research, Eindhoven, Netherlands  X. Garbet
CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France  Ph. Ghendrih  CEA, IRFM,
F-13108 Saint-Paul-lez-Durance, France  V. Grandgirard  CEA, IRFM, F-13108
Saint-Paul-lez-Durance, France  P. Manas  CEA, IRFM, F-13108 Saint-Paul-lez-
Durance, France  F. Widmer  Aix-Marseille Université, CNRS PIIM, UMR 7345
Marseille, France

(August 27, 2024)

###### Abstract

Reduced quasilinear (QL) and nonlinear (gradient-driven) models with scale
separations, commonly used to interpret experiments and to forecast turbulent
transport levels in magnetised plasmas are tested against nonlinear models
without scale separations (flux-driven). Two distinct regimes of turbulence
–either far above threshold or near marginal stability– are investigated with
Boltzmann electrons. The success of reduced models especially hinges on the
reproduction of nonlinear fluxes. Good agreement between models is found above
threshold whilst reduced models would significantly underpredict fluxes near
marginality, overlooking mesoscale flow organisation and turbulence self-
advection. Constructive prescriptions whereby to improve reduced models is
discussed.

##  I Introduction

Fusion plasmas display the property, common in dynamical systems that upon
surpassing a critical threshold, an instability may promptly build up,
inducing large fluxes which deplete the driving gradients and inhibit the
instability. Background gradients thus hover in the vicinity of nonlinear near
marginal thresholds [1]. Many strategies have been devised in modelling to
mimic natural processes. All are not equivalent and different choices may
critically affect the nature of computed statistical equilibria.

Forcing-dependent steady-states have indeed long been observed in a variety of
systems. Systems with long-range interactions, either controlled in energy or
in temperature display different equilibria [2, 3]. Swirling flows controlled
through either imposed torque or imposed velocity display distinct steady
states, as well as different dynamical regimes [4]. In magnetised fusion
plasmas, auxiliary heating and current drive force the system out of
equilibrium through a constant flux. Mimicking nature, ”flux-driven” (FD)
forcing adiabatically imparts a volumetric flux to the system whose gradients
self-consistently adapt, leading to the observation of large-scale transport
events such as avalanches [5, 6, 7, 8, 9, 10, 11], nonlinear structures such
as zonal mean flows [12, 13] or secondary patterns such as staircases [14, 15,
16, 17, 18, 19]. Such observations are absent or impaired when the system is
driven through a body force, which amounts to imposing fixed mean gradients.
This ”gradient-driven” (GD) strategy is widely used in direct numerical
computations for it is computationally efficient. It indeed enforces, with
respect to the FD approach, spatial and temporal scale separations between
equilibrium and fluctuations and solves for the fluctuations only.

Known differences between FD and GD frameworks have been documented [20, 17,
18]. Whether these are of practical incidence in fusion-relevant
configurations is non trivial. The matter is important for there are
increasing requirements for fast, reduced, yet reliable models to explore the
vast parameter space of magnetised plasma turbulence, interpret experimental
results and forecast future large experiments such as ITER. Currently, the
more advanced reduced models are based on quasilinear theory (QLT) [21, 22,
23]. With the advent of machine learning techniques, the ubiquitous closure
problem of QLT is approached through data-driven techniques that use large-
scale databases of first-principles GD computations. Systematic shortcomings,
if any, within reference GD strategies are thus likely to be carried over to
the reduced models. Comprehensive understanding of discrepancies between
first-principles FD, GD and quasilinear approaches is thus important and
timely, and the topic of this Letter.

To this end, we confront reference results from nonlinear FD gyrokinetics
using the Gysela framework [24] to state-of-the-art GD nonlinear gyrokinetics
and GD quasilinear calculations, using respectively Gkw [25] and QuaLiKiz [26,
27]. We further complement the study with confrontation with the QL flux-
driven QuaLiKiz–Jetto framework. Several instances of QuaLiKiz locally compute
at various radii flux-surface-averaged transport coefficients which are passed
on to the Jetto integrated modelling suite [28] and used to evolve profiles
through flux-driven transport equations. After a transport timescale, new
profiles are fed to QuaLiKiz and the process loops.

Main results are: (i) steady-state predictions of fluxes moderately depend on
the nature of forcing well above nonlinear threshold; (ii) near marginality
however, nonlinear and quasilinear GD models sizeably underpredict turbulent
heat transport. Under a driving flux, profiles display ’stiffness’, i.e. hover
in the vicinity of their near-threshold flux-matching values. Large, hot
devices such as ITER are expected to be stiff due to the temperature
dependence of the gyroBohm heat flux scaling, making near marginality a regime
which models must confront. Proximity to nonlinear thresholds implies
additional complexity as it favours secondary pattern formation and mesoscale
organisation. Despite this additional complexity, (iii) the underlying
assumptions of QLT hold well across nonlinear regimes. We show that transport
underprediction rather stems from the choice of closure, i.e. the nonlinear
saturation rule. This work stresses the relevance of QLT for model reductions
of turbulence whilst providing guidelines whereby reduced models can be
improved. The nonlinear and QL approaches tested here are the current
workhorse for estimating transport and confinement in turbulent fusion
plasmas. This work thus has implications for present-day experimental data
analysis and scenario extrapolation for fusion production. Novel saturation
rules should strive to incorporate near marginal flux-driven specificities,
often dubbed turbulence spreading [29, 11, 30], transport nonlocality or
staircase organisation.

##  II Two distinct regimes

Case | $\rho_{\star,50}$ | $R/a$ | $\nu_{\star,50}$ | $(q_{50}\,;\,q_{95})$ | $\tau=T_{i}/T_{e}$
---|---|---|---|---|---
near marginal | $1/250$ | $3.2$ | $0.24$ | $(1.4\,;\,4.0)$ | $1<\tau<1.3$
above thresh. | $1/350$ | $6$ | $0.02$ | $(1.7\,;\,2.8)$ | $0.9<\tau<1$
Table 1: Main plasma parameters in considered cases. Subscripts 50 and 95
respectively denote parameter values estimated at locations $r/a=0.5$ and
$r/a=0.95$.

Gysela resolves ion Larmor radius scale turbulence and collisional transport
in global tokamak geometry with no scale separation between equilibrium and
fluctuations. A centrally peaked heat source drives a deuterium plasma out of
equilibrium, which converges towards a steady temperature profile $T$ on
energy confinement times. Two paradigmatic simulation regimes are considered
in the electrostatic regime with Boltzmann electrons. Both cases are run in
the so-called “local limit” [31], at $\rho_{\ast}=\rho_{i}/a\leq 1/250$, where
comparison to local Gkw is fair. ”Local” simulations exploit to numerical
advantage the assumption of a scale separation between equilibrium and
fluctuations and solve a limited subset of the whole plasma volume which
twists around the torus due to the magnetic shear of the background magnetic
equilibrium. In this approximation, local fluxes bijectively relate to local
gradients. Here $\rho_{i}$ and $a$ are respectively the local Larmor radius
and minor radius of the torus.  The time averaged (over several tens of
turbulence correlation times in each case) radial profiles of normalised
temperature gradients $R/L_{T}=-R\partial_{r}T/T$ and zonal flow shear are
plotted in Fig.1 for the ‘Above Threshold’ (top) and ‘Near Marginal’ (bottom)
cases. The shaded areas represent temporal standard deviation. The large
deviation from the mean shearing rate in the near marginal case (bottom, right
axis) results from the meandering of staircases, which have already been
reported to play an important role in this regime of parameters [18]. The
linear (black hourglass symbols) and nonlinear (red squares) thresholds are
discussed in next section. In order to broadly span parameter regimes, main
plasma parameters vary significantly between cases, as illustrated in Table 1.

Figure 1:  Radial profiles of normalised temperature gradients (blue and red,
left axis) and zonal flow shear (green, right axis) for the ’above threshold’
(top) and ’near marginal’ (bottom) cases. Shaded areas represent temporal
standard deviation; black hourglass symbols ($\Join$) linear instability
thresholds $R/L_{T}^{\text{lin}}$ at vanishing $\mathbf{E}\times\mathbf{B}$
shear; red squares (${\color[rgb]{1,0,0}\blacksquare}$) nonlinear thresholds
$R/L_{Tc}$ estimated by Gkw in the local limit, including the Gysela
$\mathbf{E}\times\mathbf{B}$ shear.

For each set of local values of the Gysela parameters, two thresholds are to
be distinguished: linear $R/L_{T}^{\text{lin}}$ represents the normalised
temperature gradient above which an unstable mode grows, at vanishing
$\mathbf{E}\times\mathbf{B}$ shear. The inclusion of self-generated flow shear
–nonlinear in essence– introduces a second threshold $R/L_{Tc}$, indicative of
the nonlinear saturation of turbulence [32], in particular by zonal flows.
Practically, $R/L_{Tc}$ is computed in local GD frameworks as the minimum
local temperature gradient which provides a non-vanishing nonlinear heat flux.
It is estimated with Gkw through a series of increasing $R/L_{T}$ nonlinear
computations, whilst imposing local $\mathbf{E}\times\mathbf{B}$ shear values
from flux-driven Gysela. We call “above threshold” the regime such that
$R/L_{T}^{\text{lin}}<R/L_{Tc}\ll R/L_{T}$, characterised by a static, smooth
zonal mean shear $\gamma_{E}=\partial_{r}^{2}\langle\phi\rangle/B$ and
subdominant zonal fluctuations. The flux-surface average of the electric
potential is $\langle\phi\rangle$, $\tilde{\phi}$ are its fluctuations and $R$
is the tokamak major radius. We call “near marginal” regimes such that
$R/L_{T}^{\text{lin}}<R/L_{T}\lessapprox R/L_{Tc}$. Proximity from below to
$R/L_{Tc}$ is a hallmark of near marginality, featuring a well-defined
staircase pattern of flows and associated temperature corrugations which
meanders within the time interval –hence the large shear deviation– and
significant avalanching activity [18].

##  III Kubo numbers of order unity

QLT is valid [33] in the low Kubo number limit
$K=\tau_{\text{jump}}/\tau_{\text{int}}<1$, ratio of a jumping time
$\tau_{\text{jump}}$ of particles from one turbulent eddy to the next over a
nonlinear eddy-particle interaction time $\tau_{\text{int}}$. Kubo numbers are
estimated from the Gysela flux-driven data in various ways, summarised in
Tab.2. With analogy to incompressible fluids, particles trapped in turbulent
convective cells explore eddies in a typical turn-over time given by the local
vorticity $\tau_{\text{int}}^{\text{trap}}\sim
B/\langle|\nabla_{\bot}^{2}\tilde{\phi}|^{2}\rangle^{1/2}$. Whilst they
undergo this vortical motion they also drift in about
$\tau_{\text{jump}}^{\star}\sim L_{\theta}(eB/\nabla T)$ from one turbulent
structure to the next at the typical speed of the local diamagnetic velocity.
Alternatively, the slower evolution of the potential field provides a relevant
correlation time $\tau_{\text{jump}}^{{\text{corr}}}$ for turbulent
fluctuations, trade-off between unstable growth and nonlinear saturation. It
must be computed as a Lagrangian correlation time, in the co-moving frame of
the eddies to correct for the Doppler shift induced by their phase velocity.
We estimate it through image registration, following the toroidal shift
between 3-dimensional turbulence snapshots of $\tilde{\phi}$. It is compared
to turbulence-driven stochastic transport times of particles which, assuming a
diffusive ansatz for $\mathbf{E}\times\mathbf{B}$ fluctuations, drift across
eddies in about
$\tau_{\text{int}}^{{\text{diff}},x}=L_{x}/\langle|\tilde{v}_{E,x}|^{2}\rangle^{1/2}$,
with $x=\\{r,\theta\\}$ and $L_{x}$ the transverse correlation lengths
computed from Gysela outputs. Interestingly, employing a Eulerian correlation
time would result in severe Kubo number underestimation, locally up to factors
of 25.

Particle trapping $K_{\text{trap}}$ |  transverse drifts $\tau_{\text{jump}}^{\star}$
---|---
eddy turn-over $\tau_{\text{int}}^{\text{trap}}$
Random walks $K_{{\text{diff}}}^{\\{r,\theta\\}}$ |  Lagrang. correlation time $\tau_{\text{jump}}^{{\text{corr}}}$
$\mathbf{E}\times\mathbf{B}$ random walk
$\tau_{\text{int}}^{{\text{diff}},\\{r,\theta\\}}$
Table 2: Five typical wave–particle and turbulent times lead to three Kubo
number combinations, shown in Fig.2.

Figure 2:  Kubo numbers for the principal nonlinear dynamics in the problem.
Plain and plusses: turbulent radial and poloidal $\mathbf{E}\times\mathbf{B}$
velocity effect during a turbulent auto-correlation time. Circles: trapping of
particles due to turbulent vorticity during transverse crossing of the
turbulent filament.

Figure 3:  (a) and (b): gradient-driven QuaLiKiz and Gkw heat fluxes
confronted to reference flux-driven Gysela flux levels, based on the Gysela
profiles of Fig.1 and expressed in gyro-Bohm unit. Gray shaded areas represent
a standard deviation of Gysela heat fluxes during the considered time
interval, as profiles fluctuate. Red and blue shaded areas represent the
sensitivity of QuaLiKiz to these profile variations; the hourglass symbols
that of Gkw to a increased input temperature gradients. The reversed approach
is followed in panel (c): heat fluxes in QuaLiKiz–Jetto are made to match the
Gysela reference fluxes; the unknowns are the QuaLiKiz–Jetto profiles. The
remarkable agreement above threshold and large over-prediction of the
temperature gradient near marginality are consistent with results in panels
(a) and (b), emphasising reduction adequacy far from marginality and missing
physical ingredients in its vicinity.

Three Kubo numbers, combinations of the above nonlinear times are plotted in
Fig.2 for both ‘Above Threshold’ and ‘Near Marginal’ regimes. The various
definitions for $K$, coherent, provide the following picture: (i) on the basis
of order unity Kubo numbers, QLT should be marginally valid. Yet, as shown
below, key assumptions at the heart of the QLT reduction remain valid
throughout nonlinear evolution, which strengthens the case for quasi-linear
integrated modelling. Interestingly also, (ii) consistently larger $K$ values
near marginality stress the more percolative nature of transport there.
Avalanching emerges as a key theme to distinguish between ‘Above Threshold’
and ‘Near Marginal’ regimes as they likely underpin the larger $K$ values
computed near marginality. It is a likely indication that incorporating
avalanching and its zonal mean flow regulation may significantly alter
transport predictions and improve model behaviour near marginality.

##  IV Near marginal heat flux underprediction

Heat fluxes are computed with QuaLiKiz and Gkw from the Gysela time-averaged
steady-state profiles plotted in Fig.1 with the same degree of approximations
(electrostatic & Boltzmann electrons). All codes have different
normalisations, with Gysela e.g. being normalized to
$\hat{Q}^{\text{Gys}}=Q/(n_{0}T_{0}v_{Ti0})$ with
$v_{Ti0}=(T_{0}/m_{i})^{1/2}$ the ion thermal velocity. In Fig.3–(a) and (b),
fluxes from all codes are in the units of
$\hat{Q}=\hat{Q}^{\text{Gys}}/\hat{Q}^{\text{gB}}$, with
$\hat{Q}^{\text{gB}}=\rho_{\star,50}^{2}\,a/R$ and values of $\rho_{\star,50}$
and $R/a$ given in Tab.1. Consistent rescaling factors have been applied to
QuaLiKiz and Gkw fluxes. The stabilising effect of zonal flow shear is
accounted for, locally imposing the reference Gysela shear values. Non-
axisymmetric (turbulent) contributions to heat fluxes are shown in Fig.3–(a)
and (b). Turbulence spreading and profile corrugations, inherent to flux-
driven complexity are absent or hindered in QL or GD approaches. A fair
comparison thus requests that the Gysela reference data be significantly
coarse-grained before being handed over to Gkw and QuaLiKiz as input profiles
(including $\mathbf{E}\times\mathbf{B}$ shear profiles) and after
computations, when comparing e.g. fluxes. Practically, Gysela observables are
time averaged well over a typical linear growth time ($\geq
30,000\,\Omega_{ci}$), at steady-state and a radial smoothing is performed
over $60\,\rho_{i}$ to smear out visible FD specificities in the profiles —FD
turbulence correlation lengths and profile corrugations are indeed about
$10\rho_{i}$ [18]. Here, $\Omega_{ci}$ denotes the ion cyclotron frequency.
Sensitivity to gradient fluctuations, inherent to FD approaches is further
assessed by additional scans in $R/L_{T}$ and $\mathbf{E}\times\mathbf{B}$
shear within one temporal standard deviation.

Without inclusion of $\mathbf{E}\times\mathbf{B}$ shear stabilisation [i]
QuaLiKiz heat fluxes are overestimated with respect to Fig.3 by over an order
of magnitude (not shown here). With the inclusion of shear, [ii] at locations
of low or vanishing $\mathbf{E}\times\mathbf{B}$ shear and despite the
$60\,\rho_{i}$ smoothing, gradient-driven (standalone) QuaLiKiz commonly
displays [subplots (a) and (b)] variations by factors in heat fluxes from one
radial position to the next whilst Gkw exhibits much less sensitivity to
$\mathbf{E}\times\mathbf{B}$ shear stabilisation. Interestingly, this large
sensitivity of QuaLiKiz to shear is mitigated when called within the
integrated framework of Jetto [28] to allow for a flux-driven QL profile
evolution, driven by a central source that mimics the one of Gysela [subplot
(c)]. In the regime above threshold [iii] reasonable agreement in computed
fluxes is found across fidelity hierarchy. Conversely near marginality, [iv]
despite significant smoothing, heat flux discrepancies in Fig.3–(b) are well
outside allowed gradient sensitivity and fluctuation ‘error’ bars. Secular
growth of zonal flows in Gkw is responsible for the observed large flux
underprediction. This echoes previous observations casting concern on near
marginal gradient-driven predictions [17]. The soundness of separating
fluctuations from mean is thus clearly questioned near marginality.

The conclusions above are further confirmed when comparing FD Gysela to the FD
quasilinear framework of QuaLiKiz–Jetto, in the same two regimes. In subplot
(c), profiles from QuaLiKiz–Jetto are evolved until heat fluxes match the
nonlinear Gysela reference fluxes. The figure of merit now becomes how close
quasilinear profiles are at flux equilibrium with those of Gysela. Remarkable
[v] profile agreement above threshold echoes the agreement in fluxes,
displayed in panel (a). Near marginality however, [vi] a large over-prediction
of the temperature gradient is required for QuaLiKiz–Jetto to carry the same
flux as Gysela. Differently stated, both GD nonlinear fluxes and QL fluxes
(irrespective to the forcing framework) are underpredicted near marginality,
confirming that a constitutive ingredient is missing in that regime, at least
with Boltzmann electrons.

##  V Axis for improvement: saturation rules

Flux discrepancies between Gysela and Gkw likely stems from disregarding the
feedback of fluctuations on an assumed fixed ”equilibrium”. This enforces
local and single-valued flux–gradient relations, underestimates turbulence
spreading, avalanching and mesoscale organisation. All of which contribute to
transport, especially near marginality. Discrepancies with QuaLiKiz may either
come from violating assumptions central to QLT –linearity of fluctuations– or
by inheriting shortcomings akin to those of Gkw, through the choice of
saturation rules.

Figure 4:  Testing spectral distribution of key fundamental quantities: (a)
linear cross-phases, (c) heat flux, (d) saturation rule and (b) their ratio,
proxy to the dispersion relation from Gysela ($\boldsymbol{-}$), Gkw
($\times$) and QuaLiKiz ($\circ$), averaged over radial interval $0.3\leqslant
r/a\leqslant 0.7$.

To disentangle these questions, we compute in Fig.4, for all three approaches,
the complex argument $\alpha_{n}$ of the $n$-th Fourier component of the heat
flux [panel (a)] –a proxy for the linear cross-phase between electric
potential and pressure fluctuations. They depend on the toroidal wave number
$n$ labelling each eigenmode which is related to the normalised poloidal wave
vector $k_{\theta}\rho_{i}$ through:
$k_{\theta}\rho_{i}=(a/r)\,n\,q\,\rho_{\star}$. In the present cases, the
$n=25$ toroidal wave number corresponds to, at mid-radius $a/r=2$:
$k_{\theta}\rho_{i}=25*1.4*2/250=0.28$ and
$k_{\theta}\rho_{i}=25*1.7*2/350\approx 0.24$ for the Near Marginal and Above
Threshold cases, respectively. In Gysela,
$Q_{n}^{\text{Gys}}=-i|Q_{n}|\mathrm{e}^{i\alpha_{n}}=\int\hat{v}_{E,n}^{r}\frac{2}{3}\left(\hat{P}^{\ast}_{||,n}/2+\hat{P}^{\ast}_{\bot,n}\right)\mathrm{d}\theta$,
where parallel and transverse components of the pressure and the radial
component of the $\mathbf{E}\times\mathbf{B}$ drift are computed from 3D
output data. The $n$-Fourier components of the heat flux $|Q_{n}|$ and of the
squared potential $|\phi_{n}|^{2}$, i.e. the saturation rule are respectively
shown in panels (c) and (d). Panel (b) displays their ratio, a proxy for the
dispersion relation, sometimes called ‘QL flux integrals’.

Clearly, linear cross-phases display reasonable agreement, within factors of 2
and across all regimes. They are not responsible for the flux discrepancies.
Factors of disagreement –and avenues of improvement for QL modeling– are
essentially twofold: [i] $90\%$ of the heat flux is carried by modes
$n\in[\\![5,50]\\!]$ in the above threshold regime; all approaches provide
similar conclusions. Near marginality however, flux is carried in Gysela and
Gkw through $n\in[\\![5,70]\\!]$, twice the amount of active modes with
respect to QuaLiKiz. Furthermore, [ii] saturation rules are clearly
responsible for much of the observed flux discrepancies. In the above
threshold regime, flux spectra agree well [panel (c)] –though this results
from a surprising compensation: the potential, or saturation rule [panel (d)]
is under-estimated, the dispersion relation [panel (b)] is over-estimated yet
heat flux spectra in the above threshold regime are in reasonable agreement.
No such compensation occurs near marginality; there severe under-estimation of
the potential spectrum is clearly responsible for the under-prediction of
fluxes.

##  VI Outlook

The success of reduced models especially hinges on the reproduction of
nonlinear gyrokinetic fluxes [27]. Flux under-prediction in the dynamically
important regime of near marginal stability is thus a matter of importance. At
the heart of this paper lies the fact that FD and GD models provide
significantly different flux predictions with Bolzmann electrons close to
marginal stability $R/L_{Tc}$, which underpins basic discrepancies in how
nonlinear saturation of turbulence is modelled. This should foster renewed
interest in ways to complete QL or GD models near marginality. The robustness
of linear features [34] across fidelity hierarchy and across turbulent regimes
is encouraging perspective and provides constructive directions whereby
reduced models could be improved.

Discrepancies between Gysela and Gkw are strong indications that turbulence
spreading and mesoscale patterning are in fact central to accurate transport
predictions near marginality. Larger Kubo numbers near marginality, as shown
in Fig.2 reinforce this point, which is also differently stressed by recent
works in the plasma edge [11, 30]. The fact that QuaLiKiz behaves better near
marginality than Gkw is likely because the QL closure does not include the
long-lived zonal flows that quench transport near marginality in GD approaches
[17]. In QuaLiKiz, stable modes are neglected and the turbulent intensity
spectrum –chosen as a double power law in $k_{\theta}$ for unstable modes– is
fitted [35, 27] onto databases of GD nonlinear computations, similar to those
presented here with Gkw. With this procedure, QL models inherit the
shortcomings of the primitive GD models onto which they are adjusted. The near
marginal transport shortfall in QuaLiKiz indeed largely comes from issues with
the QL closure, i.e. the choice of saturation rule.

This provides directions for improvement. In physical terms, near marginal
regimes require description of transport below or at linear stability and
coupling to modes presently predicted as stable in QuaLiKiz. It also requires
to model the self-advection (spreading) of turbulent domains and the
possibility of non-monotonic flux-gradient relations. Several routes can be
explored. In the spirit of current frameworks, QL models could be trained on
near marginal flux-driven databases such as provided by the likes of Gysela.
This would likely lead to close QL models with regime-dependent turbulent
intensity spectra. Alternatively to present closures, QL models could also be
coupled to dynamic equations for the turbulence intensity, e.g. in the form of
reaction–diffusion [36] or $k-\epsilon$ equations [37], enriching accessible
nonlinear dynamics.

Ongoing studies are concerned with further characterising our “Above
Threshold” and “Near Marginal” regimes [38] when kinetic features of electron
dynamics are present. This is important to assess relevance for ITER
extrapolations. Electron dynamics is known to locally modify turbulence
organisation near low order rational surfaces [39]. Interestingly however, key
features of near marginal turbulence with Boltzmann electrons (flow
patterning, shear effectiveness and staircase organisation) which are central
here to the near marginal regime robustly endure in kinetic electron regimes
[40]. The present work provides a framework of understanding; it will be
expanded to various parameter regimes.

##  VII Acknowledgements

The authors acknowledge stimulating discussions with P.H. Diamond and
participants at the 2019 and 2021 ”Festival de Théorie” in Aix-en-Provence.
This work has been carried out within the framework of the EUROfusion
Consortium and was supported by the EUROfusion Theory and Advanced Simulation
Coordination (E-TASC) initiative under the TSVV (Theory, Simulation,
Verification and Validation) “L-H transition and pedestal physics” project
(2019–2020) and TSVV “Plasma Particle/Heat Exhaust: Gyrokinetic/Kinetic Edge
Codes” (2021–2025). It has received funding from the Euratom research and
training programme 2014-2018 and 2019-2020 under grant agreement No 633053.
The authors gratefully acknowledge funding from the European Commission
Horizon 2020 research and innovation program under Grant Agreement No. 824158
(EoCoE-II). The views and opinions expressed herein do not necessarily reflect
those of the European Commission. This research was supported in part by the
National Science Foundation under Grant No. NSF PHY-1748958. We acknowledge
PRACE for awarding us access to Joliot-Curie at GENCI@CEA, France and
MareNostrum at Barcelona Supercomputing Center (BSC), Spain.

## References

  * Diamond and Hahm [1995] P. H. Diamond and T. S. Hahm, On the dynamics of turbulent transport near marginal stability, Physics of Plasmas 2, 3640 (1995).
  * Thirring [1970] W. Thirring, Systems with negative specific heat, Zeitschrift für Physik A Hadrons and nuclei 235, 339 (1970).
  * Ellis  _et al._ [2000] R. S. Ellis, K. Haven, and B. Turkington, Large deviation principles and complete equivalence and nonequivalence results for pure and mixed ensembles, Journal of Statistical Physics 101, 999 (2000).
  * Saint-Michel  _et al._ [2013] B. Saint-Michel, B. Dubrulle, L. Marié, F. Ravelet, and F. Daviaud, Evidence for forcing-dependent steady states in a turbulent swirling flow, Phys. Rev. Lett. 111, 234502 (2013).
  * Newman  _et al._ [1996] D. E. Newman, B. A. Carreras, P. H. Diamond, and T. S. Hahm, The dynamics of marginality and self-organized criticality as a paradigm for turbulent transport, Physics of Plasmas 3, 1858 (1996).
  * Garbet and Waltz [1998] X. Garbet and R. E. Waltz, Heat flux driven ion turbulence, Physics of Plasmas 5, 2836 (1998).
  * Beyer  _et al._ [2000] P. Beyer, S. Benkadda, X. Garbet, and P. H. Diamond, Nondiffusive transport in tokamaks: Three-dimensional structure of bursts and the role of zonal flows, Physical Review Letters 85, 4892 (2000).
  * McMillan  _et al._ [2009] B. F. McMillan, S. Jolliet, T. M. Tran, L. Villard, A. Bottino, and P. Angelino, Avalanchelike bursts in global gyrokinetic simulations, Physics of Plasmas 16, 022310 (2009).
  * Idomura  _et al._ [2009] Y. Idomura, H. Urano, N. Aiba, and S. Tokuda, Study of ion turbulent transport and profile formations using global gyrokinetic full-$f$ Vlasov simulation, Nuclear Fusion 49, 065029 (2009).
  * Sarazin  _et al._ [2010] Y. Sarazin, V. Grandgirard, J. Abiteboul, S. Allfrey, X. Garbet, P. Ghendrih, G. Latu, A. Strugarek, and G. Dif-Pradalier, Large scale dynamics in flux driven gyrokinetic turbulence, Nuclear Fusion 50, 054004 (2010).
  * Singh and Diamond [2020] R. Singh and P. H. Diamond, When does turbulence spreading matter?, Physics of Plasmas 27, 042308 (2020), https://doi.org/10.1063/1.5117835 .
  * jin Kim and Diamond [2003] E. jin Kim and P. H. Diamond, Zonal flows and transient dynamics of the l-h transition, Physical Review Letters 90, 10.1103/physrevlett.90.185006 (2003).
  * Diamond  _et al._ [2005] P. H. Diamond, S.-I. Itoh, K. Itoh, and T. S. Hahm, Zonal flows in plasma—a review, Plasma Physics and Controlled Fusion 47, 35 (2005).
  * Dif-Pradalier  _et al._ [2010] G. Dif-Pradalier, P. H. Diamond, V. Grandgirard, Y. Sarazin, J. Abiteboul, X. Garbet, P. Ghendrih, A. Strugarek, S. Ku, and C. S. Chang, On the validity of the local diffusive paradigm in turbulent plasma transport, Physical Review E 82, 025401 (2010).
  * Dif-Pradalier  _et al._ [2015] G. Dif-Pradalier, G. Hornung, P. Ghendrih, Y. Sarazin, F. Clairet, L. Vermare, P. H. Diamond, J. Abiteboul, T. Cartier-Michaud, C. Ehrlacher, D. Estève, X. Garbet, V. Grandgirard, Ö. D. Gürcan, P. Hennequin, Y. Kosuga, G. Latu, P. Maget, P. Morel, C. Norscini, R. Sabot, and A. Storelli, Finding the elusive $\mathbf{E}\times\mathbf{B}$ staircase in magnetized plasmas, Physical Review Letters 114, 10.1103/physrevlett.114.085004 (2015).
  * Rath  _et al._ [2016] F. Rath, A. G. Peeters, R. Buchholz, S. R. Grosshauser, P. Migliano, A. Weikl, and D. Strintzi, Comparison of gradient and flux driven gyro-kinetic turbulent transport, Physics of Plasmas 23, 052309 (2016).
  * Peeters  _et al._ [2016] A. G. Peeters, F. Rath, R. Buchholz, Y. Camenen, J. Candy, F. J. Casson, S. R. Grosshauser, W. A. Hornsby, D. Strintzi, and A. Weikl, Gradient-driven flux-tube simulations of ion temperature gradient turbulence close to the non-linear threshold, Physics of Plasmas 23, 082517 (2016).
  * Dif-Pradalier  _et al._ [2017] G. Dif-Pradalier, G. Hornung, X. Garbet, P. Ghendrih, V. Grandgirard, G. Latu, and Y. Sarazin, The $\mathbf{E}\times\mathbf{B}$ staircase of magnetised plasmas, Nuclear Fusion 57, 066026 (2017).
  * Ashourvan  _et al._ [2019] A. Ashourvan, R. Nazikian, E. Belli, J. Candy, D. Eldon, B. A. Grierson, W. Guttenfelder, S. R. Haskey, C. Lasnier, G. R. McKee, and C. C. Petty, Formation of a high pressure staircase pedestal with suppressed edge localized modes in the diii-d tokamak, Phys. Rev. Lett. 123, 115001 (2019).
  * Nakata and Idomura [2013] M. Nakata and Y. Idomura, Plasma size and collisionality scaling of ion-temperature-gradient-driven turbulence, Nuclear Fusion 53, 113039 (2013).
  * G. Laval and Adam [2018] D. P. G. Laval and J.-C. Adam, Wave-particle and wave-wave interactions in hot plasmas: a french historical point of view, The European Physical Journal H 43, 421 (2018).
  * A.A. Vedenov and Sagdeev [1962] E. V. A.A. Vedenov and R. Sagdeev, Quasilinear theory of plasma oscillations, Nuclear Fusion Supplement 2, 465 (1962).
  * Drummond and Pines [1962] W. Drummond and D. Pines, Nonlinear stability of plasma oscillations, Nuclear Fusion Supplement 3, 1049 (1962).
  * Grandgirard  _et al._ [2016] V. Grandgirard, J. Abiteboul, J. Bigot, T. Cartier-Michaud, N. Crouseilles, G. Dif-Pradalier, C. Ehrlacher, D. Esteve, X. Garbet, P. Ghendrih, G. Latu, M. Mehrenberger, C. Norscini, C. Passeron, F. Rozar, Y. Sarazin, E. Sonnendrücker, A. Strugarek, and D. Zarzoso, A 5D gyrokinetic full-$f$ global semi-lagrangian code for flux-driven ion turbulence simulations, Computer Physics Communications 207, 35 (2016).
  * Peeters  _et al._ [2009] A. Peeters, Y. Camenen, F. Casson, W. Hornsby, A. Snodin, D. Strintzi, and G. Szepesi, The nonlinear gyro-kinetic flux tube code GKW, Computer Physics Communications 180, 2650 (2009).
  * Bourdelle  _et al._ [2007] C. Bourdelle, X. Garbet, F. Imbeaux, A. Casati, N. Dubuit, R. Guirlet, and T. Parisot, A new gyrokinetic quasilinear transport model applied to particle transport in tokamak plasmas, Physics of Plasmas 14, 112501 (2007).
  * Citrin  _et al._ [2017] J. Citrin, C. Bourdelle, F. J. Casson, C. Angioni, N. Bonanomi, Y. Camenen, X. Garbet, L. Garzotti, T. Görler, O. Gürcan, F. Koechl, F. Imbeaux, O. Linder, K. van de Plassche, P. Strand, and G. S. and, Tractable flux-driven temperature, density, and rotation profile evolution with the quasilinear gyrokinetic transport model QuaLiKiz, Plasma Physics and Controlled Fusion 59, 124005 (2017).
  * M. Romanelli [2014] e. a. M. Romanelli, Jintrac: A system of codes for integrated simulation of tokamak scenarios, Plasma and Fusion Research 9, 3403023 (2014).
  * Hahm  _et al._ [2005] T. S. Hahm, P. H. Diamond, Z. Lin, G. Rewoldt, Ö. D. Gürcan, and S. Ethier, On the dynamics of edge-core coupling, Physics of Plasmas 12, 090903 (2005).
  * Dif-Pradalier  _et al._ [2022] G. Dif-Pradalier, P. Ghendrih, Y. Sarazin, E. Caschera, F. Clairet, Y. Camenen, X. Garbet, V. Grandgirard, Y. Munschy, L. Vermare, and F. Widmer, Transport in fusion plasmas: Is the tail wagging the dog?, under review 10.21203/rs.3.rs-879691/v1 (2022).
  * Lin  _et al._ [2002] Z. Lin, S. Ethier, T. S. Hahm, and W. M. Tang, Size scaling of turbulent transport in magnetically confined plasmas, Physical Review Letters 88, 10.1103/PhysRevLett.88.195004 (2002).
  * Dimits  _et al._ [2000] A. M. Dimits, G. Bateman, M. A. Beer, B. I. Cohen, W. Dorland, G. W. Hammett, C. Kim, J. E. Kinsey, M. Kotschenreuther, A. H. Kritz, L. L. Lao, J. Mandrekas, W. M. Nevins, S. E. Parker, A. J. Redd, D. E. Shumaker, R. Sydora, and J. Weiland, Comparisons and physics basis of tokamak transport models and turbulence simulations, Physics of Plasmas 7, 969 (2000).
  * Diamond  _et al._ [2010] P. H. Diamond, S.-I. Itoh, and K. Itoh,  _Modern Plasma Physics: Volume 1, Physical Kinetics of Turbulent Plasmas_ (Cambridge University Press, 2010).
  * Besse  _et al._ [2011] N. Besse, Y. Elskens, D. F. Escande, and P. Bertrand, Validity of quasilinear theory: refutations and new numerical confirmation, Plasma Physics and Controlled Fusion 53, 025012 (2011).
  * Casati  _et al._ [2009] A. Casati, C. Bourdelle, X. Garbet, F. Imbeaux, J. Candy, F. Clairet, G. Dif-Pradalier, G. Falchetto, T. Gerbaud, V. Grandgirard, Özgur. D. Gürcan, P. Hennequin, J. Kinsey, M. Ottaviani, R. Sabot, Y. Sarazin, L. Vermare, and R. E. Waltz, Validating a quasi-linear transport model versus nonlinear simulations, Nuclear Fusion 49, 085012 (2009).
  * Heinonen and Diamond [2019] R. A. Heinonen and P. H. Diamond, Subcritical turbulence spreading and avalanche birth, Physics of Plasmas 26, 030701 (2019), https://doi.org/10.1063/1.5083176 .
  * Baschetti  _et al._ [2021] S. Baschetti, H. Bufferand, G. Ciraolo, P. Ghendrih, E. Serre, P. Tamain, and the WEST Team, Self-consistent cross-field transport model for core and edge plasma transport, Nuclear Fusion 61, 106020 (2021).
  * Garbet  _et al._ [2007] X. Garbet, Y. Sarazin, F. Imbeaux, P. Ghendrih, C. Bourdelle, O. D. Gurcan, and P. H. Diamond, Front propagation and critical gradient transport models, Physics of Plasmas 14, 122305 (2007), https://doi.org/10.1063/1.2824375 .
  * Dominski  _et al._ [2012] J. Dominski, S. Brunner, S. K. Aghdam, T. Görler, F. Jenko, and D. Told, Identifying the role of non-adiabatic passing electrons in itg/tem microturbulence by comparing fully kinetic and hybrid electron simulations, Journal of Physics: Conference Series 401, 012006 (2012).
  * Rath  _et al._ [2021] F. Rath, A. G. Peeters, and A. Weikl, Analysis of zonal flow pattern formation and the modification of staircase states by electron dynamics in gyrokinetic near marginal turbulence, Physics of Plasmas 28, 072305 (2021), https://doi.org/10.1063/5.0054358 .


11institutetext: The University of Sydney, Darlington NSW 2008, Australia
11email: {wenhua.wu, kun.hu, wenxi.yue, weiwilson.li, milena.simic,
<EMAIL_ADDRESS>22institutetext: Sydney Polytechnic Institute Pty
Ltd, Sydney, NSW, Australia
22email<EMAIL_ADDRESS>33institutetext: School of Computing,
Engineering & Mathematical Sciences, La Trobe University
33email<EMAIL_ADDRESS>
# Identity-Consistent Diffusion Network for Grading Knee Osteoarthritis
Progression in Radiographic Imaging

Wenhua Wu 11 0000-0001-7488-5487 Kun Hu Corresponding author.11
0000-0002-6891-8059 Wenxi Yue 11 0000-0002-8270-417X Wei Li 11
0000-0003-4731-3226 Milena Simic 11 Changyang Li 22 0000-0001-9400-3998 Wei
Xiang 33 0000-0001-9400-3998 Zhiyong Wang 11 0000-0002-8043-0312

###### Abstract

Knee osteoarthritis (KOA), a common form of arthritis that causes physical
disability, has become increasingly prevalent in society. Employing computer-
aided techniques to automatically assess the severity and progression of KOA
can greatly benefit KOA treatment and disease management. Particularly, the
advancement of X-ray technology in KOA demonstrates its potential for this
purpose. Yet, existing X-ray prognosis research generally yields a singular
progression severity grade, overlooking the potential visual changes for
understanding and explaining the progression outcome. Therefore, in this
study, a novel generative model is proposed, namely Identity-Consistent
Radiographic Diffusion Network (IC-RDN), for multifaceted KOA prognosis
encompassing a predicted future knee X-ray scan conditioned on the baseline
scan. Specifically, an identity prior module for the diffusion and a
downstream generation-guided progression prediction module are introduced.
Compared to conventional image-to-image generative models, identity priors
regularize and guide the diffusion to focus more on the clinical nuances of
the prognosis based on a contrastive learning strategy. The progression
prediction module utilizes both forecasted and baseline knee scans, and a more
comprehensive formulation of KOA severity progression grading is expected.
Extensive experiments on a widely used public dataset, OAI, demonstrate the
effectiveness of the proposed method.

###### Keywords:

Knee Osteoarthritis Medical Imaging Medical Assessment Generation

##  1 Introduction

Figure 1: Pipeline for KOA severity progression grading, encompassing a)
predicted future knee joint X-ray scan, and b) predicted KOA progression
severity grade.

Knee osteoarthritis (KOA) is a common form of arthritis, mainly affecting knee
joints with the symptoms of knee joint swelling, tenderness, and bone spurs.
As the major factor of physical disability, KOA has recently shown an
increasing prevalence in society, especially among the elders [9]. Currently,
there are rare treatments for curing KOA without physical surgery (total knee
replacement) [13]. Therefore, to relieve the pain and control the progression
of KOA, the diagnosis and prognosis of the disease are significant, using
different clinical assessment methods, such as medical imaging scans and
patients’ clinic questionnaires. In specific, the imaging scans are utilized
to detect the structural changes within the knee joint areas, which indicates
the characteristics of KOA. Typically, the severity of KOA is measured by the
Kellgren and Lawrence (KL) grading scale system, which encompasses five
different levels, from grade 0 (healthy cases) to grade 4 (highly severe
cases) [24].

Compared to the extensive time and human resources required for specialists to
detect changes, employing computer-aided techniques to automatically assess
the severity and progression of KOA can significantly reduce time costs. This
aids in the prompt diagnosis and prognosis of KOA at early stages [5, 20].
Methods for grading the severity and progression of KOA have evolved from
traditional machine learning approaches [35, 34] to advanced deep learning
techniques [7, 39]. These modern approaches utilize various forms of input,
including single images, multiple images, and images combined with demographic
information (such as age, gender, BMI, and the WOMAC questionnaire) [42, 27,
45, 38]. Particularly, the advancement of X-ray technology and its application
in KOA significantly improves the capability to detect and monitor the
progression [12, 18, 20]. Meanwhile, utilising the generative AI methods in
the KOA analysis field focuses on augmenting the knee radiological dataset to
improve the grading process [31]. Yet, existing X-ray prognosis research
generally yields a singular progression severity grade using the baseline
X-ray scan [18, 12, 1], overlooking the potential visual changes for
understanding and explaining the progression outcome. However, the visually
predicted outcome is also crucial in clinical practice to explain the disease
trajectory.

In addition to enhancing these methods, we hope the longitude visual content
will also aid clinicians in their decision-making by providing further
explainability and as domain knowledge to guide the algorithms. For the
methodology part, previous findings revealed that GANs yield poor bone
structures and inconsistency in patients’ identity information, which
motivated us to devise a diffusion-based method with identity information.
Identity patterns are crucial with a principle widely utilized in the forensic
field. We anticipate that incorporating such information improves the
consistency between baselines and synthesized longitude images.

Therefore, in this study, a novel deep generative model, namely Identity-
Consistent Radiographic Diffusion Network (IC-RDN), is proposed for
multifaceted KOA prognosis. IC-RDN creates a predicted future knee X-ray scan
conditioned on the baseline scan and a future KOA severity grade as
illustrated in Fig. 1. Specifically, an identity prior module for guiding the
diffusion and a downstream generation-guided progression prediction module are
introduced. Compared to a conventional image-to-image generative model,
identity priors regularize and guide the diffusion to focus more on the
clinical nuances related to the prognosis based on a contrastive learning
strategy. Using both forecasted and baseline knee scans, we anticipate a more
comprehensive and enhanced formulation of KOA severity progression measured in
grades with the progression prediction module. Extensive experiments on a
widely used public dataset, OAI [25], demonstrate the effectiveness of the
proposed method.

The key contributions of this study can be summarised in three-fold:

  * •

A novel deep generative model - IC-RDN is proposed for multifaceted KOA
prognosis prediction for the first time, encompassing radiographic visual
predictions and grading forecasting.

  * •

An identity prior module for persisting patient’s identity information on bone
X-ray images is proposed to guide the diffusion process focusing on clinical
nuances.

  * •

Comprehensive experiments with the public OAI dataset demonstrate the
effectiveness of the proposed method.

##  2 Related Works

###  2.1 Computer-Aid KOA Severity Grading

Existing studies in grading the severity of KOA can be categorised into KOA
detection and KOA multi-class grading [27, 44]. The KOA detection focuses on
categorising several or combined severity grades from the whole grading system
to achieve a coarse-grained classification. For instance, a graph neural
network has been designed to utilize the shape boundary information of knee
bones for KOA detection [42]. In addition to shape, texture features have been
independently utilized for KOA detection to examine the internal pixel
relationship [32]. To leverage both textural and shape information, the
features extracted from the shallow and deeper layers of a DenseNet model are
combined, where the shallow layers contain more shape patterns and the deeper
layers focus on textural information [27].

Multiclass grading, on the other hand, concentrates more on exploring each
grading severity level in the system. As the KL grades can be considered a
continuous measurement, with convolution neural networks, both regression and
classification strategies are studied [3, 2]. Likewise, an ordinal loss
function is proposed to formulate the continuous nature of KL grades [7]. In
parallel, the attention mechanisms identify the significant areas related to
KOA multiclass grading [11, 45]. In [22], attention mechanisms are studied to
combine and fuse multi-scale features from a network. Furthermore, to take
advantage of the symmetrical structure of knee joints, a Siamese Network is
employed to investigate the relationship between both knee joints. This
approach facilitates the analysis of intricate KOA patterns in both fully
supervised and semi-supervised manners, addressing the challenge posed by the
scarcity of annotated datasets [40, 28]. Likewise, a self-supervised method
has been introduced. This method leverages both X-ray and MRI images to expand
the pool of positive pairs in contrastive learning settings, offering a robust
solution to dataset limitations [44].

###  2.2 Computer-Aid KOA Longitude Progression Grading

Apart from KOA severity grading, predicting the progression of the disease is
also significant for medical treatment. Existing works mainly focus on
different aspects in predicting the progression of KOA: 1) Knee joint pain
status; 2) KOA risk factor prediction; 3) Joint space narrowing (JSN) and
cartilage structural change prediction; 4) KL grade prediction; and 5) total
knee replacement surgery (TKR) prediction.

Early studies on grading the progression of KOA mainly rely on conventional
machine learning algorithms. In specific, the feature selection is leveraged
in filtering the significant features contributing to the progression grades
[29]. Parallelly, the comparison among different machine learning classifiers
is conducted to select the algorithm with better progression grading
performance [43]. [8] utilises the XGBoost model to predict the knee joint
width as a KOA progression factor. To take advantage of various classification
algorithms, the ensemble strategy is deployed among RandomForest, logistic
regression, and support vector machine [23]. These machine learning-based
methods show some encouraging results in feature engineering, where the
extracted features are combined and utilized with methods such as multilayer
perceptrons (MLPs) to predict the cartilage volume from MRI scans [15].

With the powered deep learning methods explored in both general and medical
fields, end-to-end approaches have shown promising results in KOA progression
grading. An MLP is deployed as a classifier to predict progression grades [6].
Neural Networks, otherwise, are utilised to predict pain status via extracting
information from radiographics and demographics [12, 30]. Likewise, CNN-based
methods are also leveraged in analysing various KOA risk factors [6] and
predicting the JSN progression [1, 16]. To explore the KOA evolutionary
trajectory, an Adversarial Evolving Neural Network is proposed to involve the
estimated KOA patterns at different grades as auxiliary information [20].
Although predicting the direct KOA symptoms shows a clearer progression
measurement, grading the progression KL grades and predicting TKR are also
able to illustrate the progression in the clinical assessment domain. A binary
classification model is designed based on whether the patient needs to take
TKR as a future event [41, 23, 18].

Despite using medical imaging only, the information from multiple modalities,
including imaging scans, clinical questionnaires, and demographics, are
leveraged in a multimodal method to analyse the KOA progression
comprehensively [38]. However, these studies only produce a direct outcome for
progression grading ignoring the visual information from the given medical
imaging techniques. Therefore, including visual information as an additional
modality for analysing comprehensive multifaceted progression patterns is
attractive to be explored.

Figure 2: Illustration of the proposed IC-RDN method. (a) An identity prior
module formulates identity representation, which guides the diffusion process
for the generation. (b) An I2I image generative network to forecast a future
12-month X-ray scan. (c) A downstream KOA progression prediction network,
taking the baseline and the predicted 12-month X-ray scans as inputs.

##  3 Methodology

As shown in Fig. 2, the proposed IC-RDN consists of three modules: 1) an
identity prior module extracts patient’s identity representations; 2) a
diffusion-based generative module produces the future X-ray images based on
the current-stage X-ray scans and the associated identity representations; and
3) a down-stream generation-guided progression prediction module that predicts
the progression severity using both current and generated future X-ray images
to provide comprehensive multifaceted progression patterns.

###  3.1 KOA Progression Radiographic Generation

This study incorporates X-ray medical imaging scans to predict the progression
of KOA. The objective is to forecast the severity of KOA progression from two
perspectives, encompassing severity grading score and visual information
(i.e., a generated knee X-ray scan in the future - 12-month in this study).
Initially, radiographic images for joints are captured during the baseline
(first) X-ray scan, referred to as the baseline visit or 0-month visit.
Patients are then asked to undergo a follow-up X-ray scan after 12 months for
longitudinal tracking. Specifically, the $i^{\text{th}}$ longitude X-ray joint
pairs from the dataset is denoted as
$\mathbf{I}_{\text{visit}}^{i}\in\\{\mathbf{I}_{\text{baseline}}^{i},\mathbf{I}_{\text{12-month}}^{i}\\}$,
where $\mathbf{I}_{\text{visit}}^{i}\in\mathbb{R}^{C\times W\times H}$ and
$C$, $W$ and $H$ are the number of channels, the width and the height of the
X-ray joint images. For notation simplicity,
$\mathbf{I}_{\text{Baseline}}^{i}\text{ and }\mathbf{I}_{\text{12
months}}^{i}$ can be shortened as $\mathbf{I}_{\text{v00m}}^{i}\text{ and
}\mathbf{I}_{\text{v12m}}^{i}$. And the generated 12-month X-ray image for the
$i^{th}$ joint pair is denoted as ${\hat{\mathbf{I}}_{\text{12 months}}^{i}}$
or ${\hat{\mathbf{I}}_{\text{v12m}}^{i}}$. The corresponding severity grade
for the 12-month visit is indicated by $y_{KL}^{i}$ and the grades produced by
the prediction model are shown by $\hat{y}_{KL}^{i}$.

###  3.2 Diffusion-based Visual Prognosis

The diffusion-based generative module in the IC-RDN is designed to produce
12-month X-ray images from baseline visit images for understanding and
explaining the potential visual changes in terms of KOA patterns. This
incorporates a T-step Denoising Diffusion Probabilistic Model (DDPM) [19]
under an image-to-image transition manner. The DDPM consists of two primary
processes: the forward (or diffusion) process and the reverse process.
Typically, in the forward process, the original data $x_{0}\sim
q_{data}(x_{0})$ is transformed towards an isotropic Gaussian distribution.
However, within our image-to-image diffusion model, the aim is to map the data
$x_{0}\sim q_{data}(x_{0})$ to $x_{T}$, where, in our context, $x_{0}$
represents a 12-month visit X-ray image $\mathbf{I}_{\text{v12m}}^{i}$, and
$x_{T}$ corresponds to the baseline visit X-ray images
$\mathbf{I}_{\text{v00m}}^{i}$. This transformation is modelled as a fixed
Markov Chain of length $T$. Subsequently, the reverse process mirrors the
forward process, aiming to predict $x_{0}$ from $x_{T}\sim q_{data}(x_{T})$.

However, executing both processes directly in pixel space is highly resource-
intensive, necessitating substantial memory and computational power. To
mitigate this, a VQ-GAN-based encoder (see Fig. 2 (b.1.)) encodes $x_{0}$ into
a latent space $z_{0}$, wherein both the forward and reverse processes (see
Fig. 2 (b.3.)) are performed to generate $z_{T}$. A VQ-GAN-based decoder
(depicted in Fig. 2 (b.2.)) then decodes $z_{T}$ back from the latent space to
the pixel space ($x_{T}$). The objective of this stable diffusion-based
approach [33] is formulated as follows:

$\mathcal{L}_{\text{LDM}}:=\mathbb{E}_{\varepsilon(x),\varepsilon\sim\mathcal{N}(0,1),t}\left[\left\|\varepsilon-\varepsilon_{\theta}(z_{t},t)\right\|^{2}_{2}\right],$
(1)

where $\varepsilon_{\theta}(z_{t},t)$ denotes a denoising autoencoders (time-
conditional UNet) and $\theta$ contains the corresponding weights, where
$t=1...T$.

###  3.3 Knee Joint Identity Prior and Injection

To persist a subject’s identity information in generating the corresponding
12-month X-ray scan, an identity prior module is devised to extract the
identity-sensitive patterns. Specifically, we denote the identity
representation for the $i^{th}$ subject joint as
$\mathbf{F}_{\text{ID}}^{i}\in\mathbb{R}^{d}$, where $d$ indicates the
dimension of the identity representation. A CNN-based network $f(\cdot)$ is
designed to obtain the identity representation. Mathematically, we have:

$\mathbf{F}_{\text{ID}}^{i}=f(\mathbf{I}_{\text{v00m}}^{i}).$ (2)

The identity information guides and regularizes the diffusion process to not
only focus on the target visual domain distribution but also maintain
necessary information from the subjects.

For the supervision of $f$ to formulate identity information, it is proposed
to conduct training in a contrastive learning manner. A triplet loss is
introduced for this purpose as follows:

$\mathcal{L}_{\text{Triplet}}=\sum\left[\left\|f(\mathbf{I}_{\text{v00m}}^{i})-f(\mathbf{I}_{\text{v12m}}^{i})\right\|_{2}^{2}-\left\|f(\mathbf{I}_{\text{v00m}}^{i})-f(\mathbf{I}_{\text{v00m}}^{k})\right\|_{2}^{2}+\alpha\right]_{+},$
(3)

where a triplet is formulated as
$(\mathbf{I}_{\text{v00m}}^{i},\mathbf{I}_{\text{v12m}}^{i},\mathbf{I}_{\text{v00m}}^{k})$:
$\mathbf{I}_{\text{v00m}}^{i}$ indicates the anchor of current subject’s
joint; $\mathbf{I}_{\text{v12m}}^{i}$ refers to a positive sample from the
same subject’s joint; $\mathbf{I}_{\text{v00m}}^{k}$ is a negative sample with
$i\neq k$; and $\alpha$ is a margin enforcing between positive and negative
pairs. Note that the number of triplets $M$ used in the training process is
defined in the experiment implementations.

The identity information acts as a guide and regulator for the diffusion
process, ensuring that it not only aligns with the target visual domain
distribution but also preserves essential information pertaining to the
subjects. To inject the identity into the diffusion process, a cross-attention
mechanism is applied between the features from the UNet layers in
$\epsilon_{\theta}$ and $\mathbf{F}_{\text{ID}}$. For the $t^{\text{th}}$
diffusion step and $j^{\text{th}}$ layer, the attention can be formulated as
follows:

$\mathbf{A}_{j}=\text{Softmax}(\dfrac{{\mathbf{Q}}{\mathbf{K}^{\intercal}}}{\sqrt{d}})\mathbf{V},\mathbf{Q}=\phi_{j}(z_{t})\mathbf{M}^{Q}_{j}\text{,
}\mathbf{K}=\mathbf{F}_{\text{ID}}^{i}\mathbf{M}^{K}_{j}\text{, and
}\mathbf{V}=\mathbf{F}_{\text{ID}}^{i}\mathbf{M}^{V}_{j}\text{, }$ (4)

where $\mathbf{M}^{\circ_{j}}$ represents a projection matrix containing
trainable weights with a projection dimension $d$, and $\phi_{j}(z_{t})$
indicates the (flattened) intermediate representation of the UNet with
$j^{\text{th}}$ layer. The UNet can be formulated by:

$\phi_{j+1}(z_{t})=\phi_{j}(z_{t})+\mathbf{A}_{j},$ (5)

To this end, the IC-RDN can be learnt via the following scheme:

$\mathcal{L}_{\text{LDM-
Identity}}:=\mathbb{E}_{\varepsilon(x),\varepsilon\sim\mathcal{N}(0,1),t}\left[\left\|\varepsilon-\varepsilon_{\theta}(z_{t},t,f(\mathbf{I}_{\text{v00m}}^{i}))\right\|^{2}_{2}\right].$
(6)

###  3.4 Generation-Guided KOA Progression Prediction

In addition to generating 12-month X-ray image
${\hat{\mathbf{I}}_{\text{v12m}}^{i}}$ for prognosis, a singular progression
severity grade is predicted as well based on the baseline visit X-ray image
${{\mathbf{I}}_{\text{v00m}}^{i}}$ and the generated 12-month visit X-ray
image ${\hat{\mathbf{I}}_{\text{v12m}}^{i}}$. Specifically, a two-stream CNN
backbone with two branches $h_{\text{v00m}}(\cdot)$ and
$h_{\text{v12m}}(\cdot)$ is applied for the two images to produce the
associated representations:

$\mathbf{P}_{\text{v00m}}^{i}=h_{\text{v00m}}({{{\mathbf{I}}_{\text{v00m}}^{i}}}),\mathbf{P}_{\text{v12m}}^{i}=h_{\text{v12m}}({\hat{\mathbf{I}}_{\text{v12m}}^{i}}).$
(7)

The two representations are then concatenated, denoted by $\|$, and used in a
multi-layer perceptron (MLP) classifier $g(\cdot)$ to produce the predicted
progression severity grade ($\hat{y}_{KL}^{i}$), which can be formulated as:

$\hat{y}_{KL}^{i}=g(\mathbf{P}_{\text{v00m}}^{i}\|\mathbf{P}_{\text{v12m}}^{i})$
(8)

A cross-entropy loss is leveraged in the training process to guide the
prediction training process, which can be formulated as:

$\mathcal{L}_{\text{CE}}=-\sum_{i=1}^{N}y_{KL}^{i}\log\hat{y}_{KL}^{i}$ (9)

where $N$ indicates the total number of joints in the dataset.

##  4 Experimental Results and Discussion

###  4.1 Dataset and Pre-Processing

The longitudinal X-ray images used in this study were sourced from the
Osteoarthritis Initiative (OAI) [25], which is a multi-centre and longitude
study containing chronological medical imaging scans and clinical assessment
records for the subjects up to 96 months. This public dataset enables
researchers to explore the natural progression of osteoarthritis, identify
risk factors, and evaluate biomarkers for grading the KOA severity and
progression. In this study, the patients with both baseline and 12-month
visits and available 12-month severity KL grades are selected. Due to some
participants opting out of the follow-up study and missing recording of
clinical severity grading outcome, not all have recordings for the 12-month
visit image and severity KL grades. Consequently, 7,923 joints have both the
KL grade available and radiographic recordings from both the baseline and
12-month visits. As the demographic information shown in Table 1, the amount
of joints with changed severity grades among 12 months is much smaller than
that of without changes, where only 6.15% joints show a changed severity
grade. However, across these joints without the severity grade changes, the
visual structural details do not remain the same, which also shows the
importance of the proposed multifaceted progression predicting method. The
dataset was randomly divided among subjects with a distribution of 70%, 10%,
and 20% for training, validation, and testing, respectively.

Table 1: Statistics of patients, joints and KL grades in OAI longitudinal visits Longitudinal Visits | Patients | Available Joints | Joints with KL Grades
---|---|---|---
Baseline Visit | 4,551 | 9,102 | 8,953
12-Month Visit | 4,786 | 9,418 | 8,251
Longitudinal KL Changed | 487 (6.15%)
Longitudinal KL Unchanged | 7,436 (93.85%)
Both Visits | 7,923

As the raw data from the OAI contains bilateral knee joints and includes
almost all knee bones around the joint area, knee joint localisation methods
are leveraged to extract the knee joint area individually [26]. After
detecting the knee joint edge, a 224$\times$224 region surrounding the middle
of the joint was extracted from the raw data for modelling. The remaining
transformation approaches for the joint images are followed [7], including
adjusting the brightness, contrast and saturation of the images.

###  4.2 Implementation Details

A pretrained ResNet 18 on ImageNet [10] is adopted as the backbone
classification model for identity detection, which extracts the identity
features from the baseline X-ray scans. Within Identity Prior’s contrastive
learning framework, we have designated the number of triplets $M$ as 7,000 for
training and 1,000 for testing. The Identity Prior training phase is fine-
tuned using the Adam optimizer with a learning rate of 0.01, and the batch
size for this phase is established at 128. For the training process of the
diffusion-based generative module, we set the number of steps to 1,000,000
across 400 epochs. The Adam optimizer is employed here as well, with a
learning rate set at $1\times 10^{-4}$. The batch sizes for training and
testing are 64 and 8, respectively.

For the backbone models in the generation-guided KOA progression prediction
module, various CNN methods are adopted, such as VGG 16, VGG 19, ResNet 18,
ResNet 50, and DenseNet 121, for a comprehensive comparison [37, 17, 21].
Specifically, the backbone prediction models for both branches, corresponding
to $\mathbf{I}_{\text{v00m}}$ and $\hat{\mathbf{I}}_{\text{v12m}}$, share the
same architecture. In addition, the output features dimension from these
backbone models is set to 128. Then the MLP classifier $g(\cdot)$ takes the
input dimension of 256, which contains two linear layers, an activation layer,
and a dropout layer. The training process is conducted with a batch size of
16. An SGD optimizer is utilized with the learning rate of $5\times 10^{-3}$
and a weight decaying factor of $5\times 10^{-3}$ for every 5 epochs. The
entire training regimen is executed using PyTorch 2.0.1 on an NVIDIA RTX A6000
with 48 GB of video memory.

Table 2: Comparison of KOA progression grading performance for different visual generation methods. Method | IS | Progression Grading Models - Accuracy (%)
---|---|---
VGG-16 | VGG-19 | ResNet-18 | ResNet-50 | DenseNet-121
Baseline Image Only | - | 66.3 | 63.3 | 60.8 | 60.9 | 62.5
Cycle GAN [46] | 1.28 | 62.5 | 64.3 | 60.3 | 61.2 | 63.6
I2I Diffusion Model | 2.02 | 62.4 | 62.9 | 60.4 | 61.3 | 62.7
CXR-IRGen [36] | 1.76 | 61.1 | 62.7 | 60.3 | 59.1 | 61.6
PITET GAN [4] | 1.61 | 63.3 | 63.7 | 59.6 | 60.7 | 62.2
IC-RDN (ours) | 1.98 | 66.8 | 65.1 | 61.1 | 61.7 | 64.4

###  4.3 Quantitative Performance on Visual Prognosis

The evaluation of the generated future X-ray images encompasses both
quantitative and qualitative assessments. Quantitatively, we engage in
analysis through the conventional generation task metrics, inception score and
a downstream task: the prediction of KOA progression severity grade. For using
the inception score in this study, we use the Inception-v3 network trained on
the OAI 12-month visit X-ray dataset to replace the weights pre-trained on the
ImageNet [10] as the knee X-ray images are different from the natural images
at both category and structural level. As shown in Table 2, the X-ray images
generated by our proposed IC-RDN illustrate a similar IS with the X-ray images
generated by the image-to-image diffusion model and outperform the X-ray
images generated by Cycle-GAN [46] and other SOTA generative models [36, 4]
significantly. As for the similar IS between the proposed IC-RDN and the
image-to-image diffusion model, the results show similar performance on image
quality and diversity. However, the images generated by the I2I diffusion
model show a minor contribution to predicting the KOA progression severity
grades, and the challenge of identity inconsistency remains in the generated
X-ray images from the I2I diffusion model. The detailed comparison for
identity consistency is shown on Sec. 4.5.

Although the inception score is able to show the quality and diversity of the
generated future images, the downstream task of predicting the KOA progression
severity grade can be more relevant to the value of how the generated future
X-ray images contribute to exploring the KOA progression patterns. We have
introduced the approach to predict KOA severity progression (Sec. 3.4),
leveraging both the generated future X-ray image and the baseline X-ray image
to determine a progression severity grade. In addition, as stated in Sec. 4.2,
various CNN-based backbones are deployed in the generation-guided KOA
progression prediction module. The comprehensive results among these backbone
models are shown in Table. 2.

Specifically, the predicted 12-month progression X-ray images generated by
Cycle GAN, the image-to-image diffusion model, the recent SOTA generative
methods, and the proposed IC-RDN are used in the generation-guided KOA
progression grade prediction module separately. Accuracy is used as the
evaluation method to measure the performance of predicting the progression
grades. Among all backbone models used in this module, predicting the
progression severity grades using X-ray images from baseline visit and
generated by the proposed IC-RDN shows significantly outperformed results. In
detail, the image-to-image diffusion model shows a limited contribution to
providing the multifaceted KOA progression patterns. And the Cycle GAN
illustrates the minor ability to enhance the performance of predicting the KOA
progression severity grades. Therefore, the quantity evaluation demonstrates
the added value of incorporating generated future X-ray images in enhancing
the accuracy of KOA progression severity grade predictions.

Figure 3: Illustration of t-SNE analysis between IC-RDN (ours) and I2I
diffusion model.

To illustrate the performance of generated X-ray images contributing to
predicting the KOA progression grades, the t-Distributed Stochastic Neighbor
Embedding (t-SNE) technique, a sophisticated machine learning algorithm
designed for the task of high-dimensional data visualization, is utilised in
this study. In detail, to reduce the complex, multi-dimensional data into a
two-dimensional space, t-SNE converts the high-dimensional distances between
data points into conditional probabilities that represent similarities.
Specifically, given a predicted X-ray image generated by the proposed IC-RDN,
the trained VGG 19 network is used to produce a feature map corresponding to
the KOA progression grade patterns. Then, these representations are visualised
using t-SNE to show the related grades spatially. As shown in Fig. 3, the
nodes with different colours indicate the corresponding grades. In Fig. 3.a,
clear clusters for the KL grades of 0, 3, and 4 are located on the left and
right sides. The green and yellow nodes referred to KL grades 1 and 2, show
that the module for predicting the KOA progression severity grade has a minor
ability to categorise these two classes, which is the challenge among
different KOA progression prediction and severity grading methods [7, 44].
However, in Fig. 3.b, the nodes for different grades are mixed and the cluster
centres are hard to identify.

Figure 4: Comparison of generated X-ray images from Cycle GAN [46] and our
proposed IC-RDN on knee joints i-viii.  Figure 5: Generative results from IC-
RDN among different patients’ joints from i to xxii. The detailed joint
structure comparison between baseline visit images, ground truth and results
from IC-RDN shows our model’s ability to generate progression patterns.

###  4.4 Qualitative Performance on Visual Prognosis

In terms of qualitative analysis, the comparison between X-ray images
synthesized by the Cycle GAN and those produced by the proposed IC-RDN offers
a compelling visual assessment, as depicted in Fig. 4. The comparative study
highlights notable distortions in bone structures emanating from the Cycle GAN
approach, particularly in the articulation of joint spaces. In stark contrast,
our proposed IC-RDN exhibits a pronounced fidelity in preserving the integrity
of key bone edges and the intricacy of inner textures, thereby underscoring
the superior capability of IC-RDN in maintaining structural information.

Further elucidation is provided through an expanded corpus of X-ray images
generated by the IC-RDN, which is shown in Fig. 5. This compilation serves not
only to reinforce the initial observations regarding the preservation of
critical structural details but also to demonstrate the consistent performance
of the IC-RDN across a diverse set of examples. The consistent replication of
complex anatomical features by the IC-RDN model is indicative of its
sophisticated understanding of underlying bone structures, which is paramount
for illustrating the multifaceted KOA progression patterns.

Moreover, this comparison elucidates the inherent limitations of traditional
GAN-based approaches in contexts where the preservation of exact structural
features is paramount. The propensity for distortion within Cycle GAN-
generated images can significantly impact the clinical utility of such
synthesized X-ray images, particularly in monitoring the KOA progression. In
contrast, the proposed IC-RDN shows its adeptness at capturing and reproducing
the nuanced details of bone structure suggesting its potential as a
multifaceted KOA progression prediction outcome.

Figure 6: Examples of generated X-ray images with identity consistency for
joint i-x. The first line (outputs from the Image-to-Image Diffusion Model)
shows the generative results without Identity Prior. The second line indicates
the ground truth images. The third line refers to the results generated by our
IC-RDN approach.

###  4.5 Analysis on Identity Consistency

To further evaluate the efficacy of the identity prior module, we delve into a
comprehensive comparison involving the X-ray images generated by the I2I
diffusion model and the proposed IC-RDN. The visual examples presented in Fig.
6 meticulously outline the distinctions between these methods alongside the
actual 12-month visit X-ray images, where the key differences are highlighted
by the bounding boxes. Initially, we observe the generated X-ray images devoid
of identity persistency, which starkly contrasts with the subsequent imagery.
A noteworthy observation is the preservation of identity features in the
X-rays produced by our IC-RDN, as evidenced by the distinct texture and
precision of bone edge structures. This preservation extends to intricate
details within the knee joint, such as the orientation of the menisci and the
nuanced protrusions along the bone edges, which are critical for maintaining
identity persistency.

###  4.6 Ablation Study

To better evaluate the contribution of the proposed Identity Prior module, an
ablation study is conducted by removing the Identity Prior and only the
Diffusion backbone remaining. The results in Table 3 demonstrate the
effectiveness of the Identity Prior module with a superior longitude
prediction accuracy among all selected classification methods.

Table 3: Ablation study of the proposed IC-RDN Method | IS | Progression Grading Models - Accuracy (%)
---|---|---
VGG16 | VGG19 | ResNet18 | ResNet50 | DenseNet121
IC-RDN w/o Identity Prior | 1.72 | 61.6 | 61.7 | 60.2 | 61.0 | 62.5
IC-RDN (ours) | 1.98 | 66.8 | 65.1 | 61.1 | 61.7 | 64.4

###  4.7 Limitations

The major limitation of this study regards the generation quality. Directly
using X-ray scans generated by IC-RDN for a KOA severity grading showcases a
suboptimal result compared with using the baseline visit X-ray to predict the
KOA progression grade directly. Though the predicted 12-month scan can enhance
the KOA severity prediction with using the baseline visit jointly. In
addition, computer hallucinations [14] are a crucial issue among large model
fields. Particularly, the analysis regarding the disease-specific changes is
insufficient. Therefore, careful consideration of applying the generative
method to actual clinical usage is needed.

##  5 Conclusion

In this study, a generative-based method, IC-RDN is proposed for KOA
prognosis, which forecasts future X-ray scans for joints and singular
progression KOA severity grades. This obtains multifaceted KOA progression
patterns for understanding and explaining KOA prognosis. Specifically, an
identity prior module for the diffusion and a downstream generation-guided
progression prediction module are introduced for regularising and guiding the
diffusion to focus more on the clinical nuances related to the prognosis and
producing a more precise KOA severity progression grading. Extensive
experiments on a widely used public dataset, OAI, demonstrate the
effectiveness of the proposed method. In our future studies, the
disentanglement of the representations related to the identity and KOA
severity will be explored for fine-grained guidance in the X-ray image
generation process. Meanwhile, the multimodal patterns among different medical
imaging scans (_e.g_. X-ray, MRI, and CT) and patients’ demographic
information will be involved in producing a fine-grained KOA progression
assessment outcome.

## Acknowledgements

This study was partially supported by Australian Research Council (ARC) grant
DP210102674. The OAI [25] is a public-private partnership comprised of five
contracts (N01-AR-2-2258; N01-AR2-2259; N01-AR-2- 2260; N01-AR-2-2261;
N01-AR-2-2262) funded by the National Institutes of Health, a branch of the
Department of Health and Human Services, and conducted by the OAI Study
Investigators.

## References

  * [1] Almhdie-Imjabbar, A., Nguyen, K.L., Toumi, H., Jennane, R., Lespessailles, E.: Prediction of knee osteoarthritis progression using radiological descriptors obtained from bone texture analysis and siamese neural networks: data from OAI and most cohorts. Arthritis Research & Therapy 24(1), 1–13 (2022)
  * [2] Antony, J., McGuinness, K., Moran, K., O’Connor, N.E.: Automatic detection of knee joints and quantification of knee osteoarthritis severity using convolutional neural networks. In: International Conference on Machine Learning and Data Mining. vol. 10358, pp. 376–390. Springer (2017)
  * [3] Antony, J., McGuinness, K., O’Connor, N.E., Moran, K.: Quantifying radiographic knee osteoarthritis severity using deep convolutional neural networks. In: Int. Conf. Pattern Recog. pp. 1195–1200. IEEE (2016)
  * [4] Berman, O., Oz, N., Mendlovic, D., Sochen, N., Cohen, Y., Klapp, I.: PETIT-GAN: physically enhanced thermal image-translating generative adversarial network. In: 2024 IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 1607–1616 (2024)
  * [5] Brahim, A., Jennane, R., Riad, R., Janvier, T., Khedher, L., Toumi, H., Lespessailles, E.: A decision support tool for early detection of knee osteoarthritis using x-ray imaging and machine learning: data from the Osteoarthritis Initiative. Computerized Medical Imaging and Graphics 73, 11–18 (2019)
  * [6] Chan, L., Li, H., Chan, P., Wen, C.: A machine learning-based approach to decipher multi-etiology of knee osteoarthritis onset and deterioration. Osteoarthritis and Cartilage Open 3(1), 100135 (2021)
  * [7] Chen, P., Gao, L., Shi, X., Allen, K., Yang, L.: Fully automatic knee osteoarthritis severity grading using deep neural networks with a novel ordinal loss. Comput. Med. Imag. Grap. 75, 84–92 (2019)
  * [8] Cheung, J.C.W., Tam, A.Y.C., Chan, L.C., Chan, P.K., Wen, C.: Superiority of multiple-joint space width over minimum-joint space width approach in the machine learning for radiographic severity and knee osteoarthritis progression. Biology 10(11), 1107 (2021)
  * [9] Cui, A., Li, H., Wang, D., Zhong, J., Chen, Y., Lu, H.: Global, regional prevalence, incidence and risk factors of knee osteoarthritis in population-based studies. EClinicalMedicine 29 (2020)
  * [10] Deng, J., Dong, W., Socher, R., Li, L.J., Li, K., Fei-Fei, L.: Imagenet: a large-scale hierarchical image database. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 248–255. IEEE (2009)
  * [11] Górriz, M., Antony, J., McGuinness, K., Giró-i Nieto, X., O’Connor, N.E.: Assessing knee OA severity with CNN attention-based end-to-end architectures. In: International Conference on Medical Imaging with Deep Learning. pp. 197–214. PMLR (2019)
  * [12] Guan, B., Liu, F., Mizaian, A.H., Demehri, S., Samsonov, A., Guermazi, A., Kijowski, R.: Deep learning approach to predict pain progression in knee osteoarthritis. Skeletal Radiology pp. 1–11 (2022)
  * [13] Guermazi, A., Hunter, D.J., Roemer, F.W.: Plain radiography and magnetic resonance imaging diagnostics in osteoarthritis: validated staging and scoring. J. Bone Joint Surg. Am. 91(Supplement_1), 54–62 (2009)
  * [14] Gunjal, A., Yin, J., Bas, E.: Detecting and preventing hallucinations in large vision language models. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 38, pp. 18135–18143 (2024)
  * [15] Hafezi-Nejad, N., Guermazi, A., Roemer, F.W., Hunter, D.J., Dam, E.B., Zikria, B., Kwoh, C.K., Demehri, S.: Prediction of medial tibiofemoral compartment joint space loss progression using volumetric cartilage measurements: data from the FNIH OA biomarkers consortium. European Radiology 27, 464–473 (2017)
  * [16] Halilaj, E., Le, Y., Hicks, J.L., Hastie, T.J., Delp, S.L.: Modeling and predicting osteoarthritis progression: data from the Osteoarthritis Initiative. Osteoarthritis and Cartilage 26(12), 1643–1650 (2018)
  * [17] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: IEEE Conference on Computer Vision and Pattern Recognition. pp. 770–778 (2016)
  * [18] Hirvasniemi, J., Runhaar, J., van der Heijden, R., Zokaeinikoo, M., Yang, M., Li, X., Tan, J., Rajamohan, H., Zhou, Y., Deniz, C., et al.: The knee osteoarthritis prediction (KNOAP2020) challenge: An image analysis challenge to predict incident symptomatic radiographic knee osteoarthritis from MRI and x-ray images. Osteoarthritis and Cartilage 31(1), 115–125 (2023)
  * [19] Ho, J., Jain, A., Abbeel, P.: Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems 33, 6840–6851 (2020)
  * [20] Hu, K., Wu, W., Li, W., Simic, M., Zomaya, A., Wang, Z.: Adversarial evolving neural network for longitudinal knee osteoarthritis prediction. IEEE Transactions on Medical Imaging 41(11), 3207–3217 (2022)
  * [21] Huang, G., Liu, Z., Van Der Maaten, L., Weinberger, K.Q.: Densely connected convolutional networks. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 4700–4708 (2017)
  * [22] Jain, R.K., Sharma, P.K., Gaj, S., Sur, A., Ghosh, P.: Knee osteoarthritis severity prediction using an attentive multi-scale deep convolutional neural network. Multimedia Tools and Applications 83(3), 6925–6942 (2024)
  * [23] Jamshidi, A., Pelletier, J.P., Labbe, A., Abram, F., Martel-Pelletier, J., Droit, A.: Machine learning–based individualized survival prediction model for total knee replacement in osteoarthritis: data from the Osteoarthritis Initiative. Arthritis Care & Research 73(10), 1518–1527 (2021)
  * [24] Kellgren, J.H., Lawrence, J.: Radiological assessment of osteo-arthrosis. Annals of the Rheumatic Diseases 16(4), 494 (1957)
  * [25] Lester, G.: Clinical research in OA-the NIH Osteoarthritis Initiative. J. Musculoskelet Neuronal Interact 8(4), 313–314 (2008)
  * [26] Lindner, C., Thiagarajah, S., Wilkinson, J.: Fully automatic segmentation of the proximal femur using random forest regression voting. IEEE Trans. Med. Imaging 32, 1462–1472 (2013)
  * [27] Nasser, Y., Hassouni, M.E., Jennane, R.: Discriminative deep neural network for predicting knee osteoarthritis in early stage. In: International Workshop on PRedictive Intelligence In MEdicine. pp. 126–136. Springer (2022)
  * [28] Nguyen, H.H., Saarakkala, S., Blaschko, M.B., Tiulpin, A.: Semixup: in-and out-of-manifold regularization for deep semi-supervised knee osteoarthritis severity grading from plain radiographs. IEEE Trans. Med. Imaging 39(12), 4346–4356 (2020)
  * [29] Ntakolia, C., Kokkotis, C., Moustakidis, S., Tsaopoulos, D.: Prediction of joint space narrowing progression in knee osteoarthritis patients. Diagnostics 11(2), 285 (2021)
  * [30] Pierson, E., Cutler, D.M., Leskovec, J., Mullainathan, S., Obermeyer, Z.: An algorithmic approach to reducing unexplained pain disparities in underserved populations. Nature Medicine 27(1), 136–140 (2021)
  * [31] Prezja, F., Annala, L., Kiiskinen, S., Ojala, T.: Exploring the efficacy of base data augmentation methods in deep learning-based radiograph classification of knee joint osteoarthritis. Algorithms 17(1), 8 (2023)
  * [32] Ribas, L.C., Riad, R., Jennane, R., Bruno, O.M.: A complex network based approach for knee osteoarthritis detection: data from the Osteoarthritis Initiative. Biomedical Signal Processing and Control 71, 103133 (2022)
  * [33] Rombach, R., Blattmann, A., Lorenz, D., Esser, P., Ommer, B.: High-resolution image synthesis with latent diffusion models. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 10684–10695 (2022)
  * [34] Shamir, L., Ling, S.M., Scott, W., Hochberg, M., Ferrucci, L., Goldberg, I.G.: Early detection of radiographic knee osteoarthritis using computer-aided analysis. Osteoarthritis and Cartilage 17(10), 1307–1312 (2009)
  * [35] Shamir, L., Ling, S.M., Scott, W.W., Bos, A., Orlov, N., Macura, T.J., Eckley, D.M., Ferrucci, L., Goldberg, I.G.: Knee x-ray image analysis method for automated detection of osteoarthritis. IEEE Transactions on Biomedical Engineering 56(2), 407–415 (2008)
  * [36] Shentu, J., Al Moubayed, N.: Cxr-irgen: an integrated vision and language model for the generation of clinically accurate chest x-ray image-report pairs. In: 2024 IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 5200–5209 (2024)
  * [37] Simonyan, K., Zisserman, A.: Very deep convolutional networks for large-scale image recognition. In: International Conference on Learning Representations (2015)
  * [38] Tiulpin, A., Klein, S., Bierma-Zeinstra, S.M., Thevenot, J., Rahtu, E., Meurs, J.v., Oei, E.H., Saarakkala, S.: Multimodal machine learning-based knee osteoarthritis progression prediction from plain radiographs and clinical data. Scientific Reports 9(1), 20038 (2019)
  * [39] Tiulpin, A., Saarakkala, S.: Automatic grading of individual knee osteoarthritis features in plain radiographs using deep convolutional neural networks. Diagnostics 10(11), 932 (2020)
  * [40] Tiulpin, A., Thevenot, J., Rahtu, E., Lehenkari, P., Saarakkala, S.: Automatic knee osteoarthritis diagnosis from plain radiographs: a deep learning-based approach. Sci. Rep. 8(1), 1–10 (2018)
  * [41] Tolpadi, A.A., Lee, J.J., Pedoia, V., Majumdar, S.: Deep learning predicts total knee replacement from magnetic resonance images. Scientific Reports 10(1), 6371 (2020)
  * [42] von Tycowicz, C.: Towards shape-based knee osteoarthritis classification using graph convolutional networks. In: IEEE International Symposium on Biomedical Imaging. pp. 750–753. IEEE (2020)
  * [43] Widera, P., Welsing, P.M., Ladel, C., Loughlin, J., Lafeber, F.P., Petit Dop, F., Larkin, J., Weinans, H., Mobasheri, A., Bacardit, J.: Multi-classifier prediction of knee osteoarthritis progression from incomplete imbalanced longitudinal data. Scientific Reports 10(1), 8427 (2020)
  * [44] Wu, W., Hu, K., Yue, W., Li, W., Simic, M., Li, C., Xiang, W., Wang, Z.: Self-supervised multimodal fusion network for knee osteoarthritis severity grading. In: International Conference on Digital Image Computing: Techniques and Applications. pp. 57–64. IEEE (2023)
  * [45] Zhang, B., Tan, J., Cho, K., Chang, G., Deniz, C.M.: Attention-based CNN for KL grade classification: data from the Osteoarthritis Initiative. In: International Symposium on Biomedical Imaging. pp. 731–735. IEEE (2020)
  * [46] Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired image-to-image translation using cycle-consistent adversarial networks. In: Proceedings of the IEEE International Conference on Computer Vision. pp. 2223–2232 (2017)