Stochastic gradient Markov Chain Monte Carlo algorithms are popular samplers for approximate inference, but they are generally biased. We show that many recent versions of these methods (e.g. [Chen et al., 2014]) cannot be corrected using Metropolis-Hastings rejection sampling, because their acceptance probability is always zero. We can fix this by employing a sampler with realizable backwards trajectories, such as Gradient-Guided Monte Carlo [Horowitz, 1991], which generalizes stochastic gradient Langevin dynamics [Welling and Teh, 2011] and Hamiltonian Monte Carlo. We show that this sampler can be used with stochastic gradients, yielding nonzero acceptance probabilities, which can be computed even across multiple steps.

§ INTRODUCTION

Exact posterior inference in many Bayesian models is intractable, which has led to the development of different approximate inference schemes.
One popular class of such algorithms are Markov Chain Monte Carlo (MCMC) methods. These simulate a Markov chain, constructed so that its stationary distribution is the true posterior [Geyer, 1992].
Among the MCMC methods, Hamiltonian Monte Carlo (HMC) has achieved recognition as being the “gold standard”, mostly for its favorable convergence properties [Neal, 2012].

While HMC yields great performance, its adoption in many real applications is hindered by the fact that it requires passes over the whole data set to compute one gradient step, which is prohibitively expensive in many settings [Chen et al., 2014].
Especially for large-data Bayesian models such as Bayesian neural networks (BNNs), stochastic gradient MCMC (SG-MCMC) methods are more popular, due to their favorable computational costs and their similarity to stochastic gradient descent (SGD).
Arguably, the most popular SG-MCMC method for BNNs in recent years has been stochastic gradient HMC <cit.>, although its predecessor stochastic gradient Langevin dynamics <cit.> also enjoys considerable following.

Since all of these approaches are approximate, in order to ensure that they sample from the true posterior, rejection sampling is often necessary, for instance, following the Metropolis-Hastings (MH) scheme [Metropolis et al., 1953, Hastings, 1970].
However, in this work, we draw attention to the unfortunate fact that SGHMC does not allow for MH rejection sampling, since its acceptance probability is always zero.
We then revisit gradient-guided Monte Carlo (GGMC) [Horowitz, 1991], which can overcome this problem.
We show that GGMC provides a unified framework for SGLD and SGHMC and that it yields nonzero acceptance probabilities, which can even be computed across multiple steps.

Our contributions are as follows:

* We show that the SG-MCMC schemes used in recent work [Chen et al., 2014, Wenzel et al., 2020] always have acceptance probability zero, because the backward trajectory is not realizable in the symplectic Euler-Maruyama integrator they employ (Sec. <ref>).

* We revisit GGMC, a variant of HMC. It generalises HMC and SGLD, always has positive acceptance probability and, with correction, accepts stochastic gradients (Sec. <ref>). The scheme was independently discovered in the statistics [Neal, 2012, Horowitz, 1991] and molecular dynamics literatures [Bussi and Parrinello, 2007, Leimkuhler and Matthews, 2015]. The latter have shown that it simulates Langevin dynamics.

* We examine an augmentation scheme to calculate the acceptance probabilities across multiple stochastic steps (<ref>, <ref>; [Zhang et al., 2020]). This allows us to use stochastic gradients, but sample from the exact posterior, by doing one MH correction using the exact likelihood after many stochastic gradient steps.

Based on our findings, we believe that work in Bayesian deep learning should discard the Euler-Maruyama scheme used by SGHMC and its variants [Ding et al., 2014, Ma et al., 2019]. Instead, to have confidence in our sampling, we should use reversible integrators <cit.> and keep track of their acceptance probability.

§ BACKGROUND: HAMILTONIAN MONTE CARLO

We wish to sample from a target distribution $\pi\bra{\vtheta} \propto \Ptarg{\vtheta}$, which we can
only evaluate up to a normalizing constant. HMC first augments the parameter
space $\vtheta \in \mathbb{R}^d$ with a momentum $\vm\in\mathbb{R}^{d}$ of the
same dimension. It constructs a Markov chain that alternates between re-sampling the momentum and then
following an ordinary differential equation (ODE), specified by the Hamiltonian $H(\vtheta, \vm)$,
for a given amount of time.
The stationary distribution of this Markov chain (with exact simulation) is the Boltzmann distribution

\begin{equation}
\Pboltzmann{\vtheta, \vm} \propto \exp\bra{-\frac{1}{T}H(\vtheta, \vm)} = \exp\bra{-\frac{1}{T} U(\vtheta)} \exp\bra{-\frac{1}{T}K(\vm)}.
\label{eq:boltzmann}
\end{equation}

Here, $U(\vtheta)$ is the potential energy, which we relate to the unnormalised target distribution by $U(\vtheta) = -\log \Ptarg{\vtheta}$.
The kinetic energy $K(\vm) = \frac{1}{2}\vm^{\tp}\mM^{{-1}}\vm$  makes the momenta $\vm$ Gaussian, and thus easy to sample.
$\vM$ is a positive
definite (usually diagonal) preconditioner, the mass matrix.
$T \in \mathbb{R}^{+}$ is a positive temperature. For HMC, $T=1$, but we will relax this constraint in the next section. By construction, the tempered target $\Ptarg{\vtheta}^{1/T}$ is a marginal of the Boltzmann distribution.
We may write the ODE's full state as $\vs = (\vtheta, \vm)$.

Unfortunately, we cannot exactly simulate the ODE for most posteriors, so we must discretise it with an integrator. This unavoidably introduces errors into the simulation, which have to be controlled if we want to guarantee that the sampling is accurate.

The standard way to control these errors is via the Metropolis-Hastings acceptance probability. For this, we define the transition kernel: $\Ptrans{\vs_*}{\vs_n} = 1$ if the approximate simulation starting at $\vs_n=(\vtheta_n, \vm_n)$ ends at $\vs_*=(\vtheta_*, \vm_*)$, otherwise 0.

Then, we accept the new sample (i.e., set $\vs_{n+1}=\vs_* $) with probability $\min\bra{1, \mhaccept{\vs_{*}}{\vs_{n}}}$, where

\begin{equation}
\mhaccept{\vs_{*}}{\vs_{n}} = \frac{\Pboltzmann{\vtheta_*, \vm_*}}{\Pboltzmann{\vtheta_n, \vm_n}} \frac{\Ptrans{\vtheta_n, -\vm_n}{\vtheta_*, -\vm_*}}{\Ptrans{\vtheta_*, \vm_*}{\vtheta_n, \vm_n}},
\label{eq:metropolis-hastings}
\end{equation}

otherwise we reject the sample, that is, set $\vs_{n+1}=\vs_{n}$. We negate the momentum so the integrator applied backwards would retrace the same steps.[Hamiltonian dynamics are time-reversible. We also need to use a time-reversible integrator, like the leapfrog or velocity Verlet, so the ratio of $g$ in <ref> is 1.]
Afterwards, we are free to negate the momentum again without a rejection step, since $K(-\vm) = K(\vm)$, so it leaves the target distribution (<ref>) invariant <cit.>. For regular HMC, this can be ignored because the momentum gets resampled immediately.

§ SGHMC ACCEPTANCE PROBABILITIES ARE ZERO

All the Markov chains we consider in this paper are discretisations of
Langevin's stochastic differential equations (SDEs) which, borrowing notation
from Wenzel et al., 2020, are

\begin{align}
\dd\vtheta &= \vM^{{-1}} \vm \,\dd{t}
\label{eq:langevin-parameter} \\
\dd\vm &= -\gamma \vm \,\dd{t} - \gradat{}\,\dd{t} + \vM^{1/2}\sqrt{2\gamma T}\,\dd{\vW} .
\label{eq:langevin-momentum}
\end{align}

Here, $\gamma \in \mathbb{R}^{+}$ is a friction coefficient and $\vW$ is a standard
Wiener process, with covariance $\mI$.

If $\gamma = 0$, these become Hamilton's equations, used for HMC.
Setting $\dd\vm=0$, we obtain the overdamped Langevin equation, which discretises
to the Metropolis-adjusted Langevin algorithm (MALA; e.g., [Roberts and Rosenthal, 1998]).

For time step $h$, the symplectic Euler-Maruyama scheme in <cit.> and Wenzel et al., 2020 simulates these SDEs as follows. Using noise $\vepsilon \sim \Normal{{\boldsymbol 0}, \mI}$,

\begin{align}
\vm_{n+1} &= (1 - h\gamma) \vm_{n} - h\gradat{_{n}} + \sqrt{h}\mM^{1/2}\sqrt{2\gamma T} \vepsilon
\label{eq:symplectic-euler-momentum} \\
\vtheta_{n+1} &= \vtheta_{n} + h\mM^{{-1}} \vm_{n+1} \;.
\label{eq:symplectic-euler-parameter}
\end{align}

Like in HMC, it is the discretisation, not the SDE, that determines our simulation's Markov transition kernel $g$. Unfortunately, we will see that the backward transition is unrealizable, so its density $\Ptrans{\vtheta_{n}, -\vm_{n}}{\vtheta_{n+1}, -\vm_{n+1}}$ is almost always zero. Because this density appears in the numerator of <ref>, the resulting acceptance probability for the Euler-Maruyama scheme is zero as well.

For any starting point $(\vtheta_{n},\vm_{n})$ and step size $h > 0$, we sample
$(\vtheta_{n+1},\vm_{n+1})$ using
equations (<ref>, <ref>).
Then, the backward transition density
$\Ptrans{\vtheta_{n}, -\vm_{n}}{\vtheta_{n+1}, -\vm_{n+1}}$ is zero, with probability 1 over the randomness of the forward transition.

Consider the transition distribution $\Ptrans{\vtheta_*, \vm_*}{\vtheta, \vm}$. Because it has been generated from equations (<ref>, <ref>), any realization from $g$ must satisfy them. In particular, it has to satisfy <ref>, that is, $\vtheta_{*} = \vtheta + h\mM^{-1}\vm_{*}$. Thus,

for the purposes of the acceptance probability (eq. <ref>), we can define the density of $g$ as zero outside of its support,

\begin{equation}
\Ptrans{\vtheta_{*}, \vm_{*}}{\vtheta,\vm} = \begin{cases}
\Normal{\vm_{*} \mvbar (1 - h\gamma)\vm - h\gradat{}, h\mM2\gamma T}
&\text{if } \vtheta_{*} = \vtheta + h\mM^{-1}\vm_{*} \\
0 & \text{otherwise}.\end{cases}
\label{eq:transition-density}
\end{equation}

Now, let us check whether this condition can simultaneously hold for the  forward and backward transitions in <ref>. That is, is it possible that

\begin{equation}
\vtheta_{n+1} = \vtheta_n + h\mM^{-1}\vm_{n+1}
\quad\text{and}\quad
\vtheta_{n} = \vtheta_{n+1} + h\mM^{-1}(-\vm_{n}) \; \text{?}
\label{eq:forward-and-backward-conditions}
\end{equation}

Solving the linear system we can see that, since $\mM$ is positive definite
and $h>0$, <ref> is true if and only if
$\vm_n = \vm_{n+1}$. Since the step size is $h > 0$, this only happens when
the Gaussian random draw $\vepsilon$ in <ref> is
equal to
$\vv = \bra{h\mM 2\gamma T}^{-1/2} \bra{ h\gamma \vm_n + h\gradat{_n}}$.
Because the mean of $\vepsilon$ is zero (and thus $\ne \vv$) and its
distribution is Gaussian, this is a probability zero event.[This also holds for any non-singular distributions on $\vepsilon$. It is also undesirable, because the momentum is supposed to change over time when following Langevin dynamics.]

Thus, with probability 1, $\vepsilon \ne \vv$, which implies $\vm_n \ne \vm_{n+1}$. We know the left-hand-side expression of <ref> is true, so that implies $\vtheta_n \ne \vtheta_{n+1} + h\mM^{-1}(-\vm_n)$, which implies by <ref> that $\Ptrans{\vtheta_{n}, -\vm_{n}}{\vtheta_{n+1}, -\vm_{n+1}} = 0$.

The easiest way to ensure $\vm_{n} = \vm_{n+1}$, such that the acceptance probability
remains nonzero, is to set the time step $h=0$. This, however, makes the sampler
useless. If we take the limit $h \to 0$ instead, the Euler-Maruyama scheme gets
arbitrarily close to the true SDE trajectory. Crucially however, per <ref>, the acceptance probability
remains $0$ for any $h>0$. Thus, it is impossible to use the acceptance probability to monitor the discretisation error. The Euler-Maruyama scheme cannot satisfy detailed balance. So far, we have not considered stochastic gradients, but this result includes them: it still holds if we substitute an arbitrary $\vg_{n}$ for $\gradat{_n}$.

§ GGMC SOLVES THE PROBLEM BY USING A DIFFERENT INTEGRATOR

To get a nonzero acceptance probability, we can just use a backward-realizable integrator. The one we present here has been discovered several times in statistics [Horowitz, 1991, Neal, 1993] and molecular dynamics [Bussi and Parrinello, 2007, Leimkuhler and Matthews, 2015]. We call it Gradient-guided Monte Carlo (GGMC), following Horowitz, 1991, because it operates similarly to random walk Metropolis-Hastings on the momentum, but with a parameter gradient influencing the random walk's direction.
<cit.> calls it the “stochastic dynamics” method, or “Langevin Monte Carlo (LMC) with partial momentum refreshment” [Neal, 2012].
Leimkuhler and Matthews, 2015 call this integrator “OBABO” because of the order of its steps.

The integrator splits Langevin dynamics
(eqs. <ref>–<ref>) into three parts,
which can be solved exactly, and applies them in sequence [Leimkuhler and Matthews, 2015]. Let $a = e^{-\gamma h}$, and $\vepsilon,\vepsilon' \sim \Normal{{\boldsymbol 0}, \mI}$ be Gaussian random variables with identity covariance. Then

\begin{align}
\tag{O.1} \label{eq:O.1}
\vm_{n+{1/4}} &= \sqrt{a} \vm_{n} + \sqrt{(1 - a) T} \mM^{1/2}\vepsilon, \\
\tag{B.1}\label{eq:B.1}
\vm_{n+{1/2}} &= \vm_{n+{1/4}} - \frac{h}{2}\gradat{_{n}}, \\
\tag{A}\label{eq:A}
\vtheta_{n+1} &= \vtheta_{n} + h\mM^{{-1}}\vm_{n+{1/2}}, \\
\tag{B.2}\label{eq:B.2}
\vm_{n+{3/4}} &= \vm_{n+{1/2}} - \frac{h}{2}\gradat{_{n+1}}, \\
\tag{O.2}\label{eq:O.2}
\vm_{n+1} &= \sqrt{a} \vm_{n+{3/4}} + \sqrt{(1 - a) T} \mM^{{1/2}} \vepsilon'.
\end{align}

Steps <ref> and <ref> are partial momentum refreshment, which is exact; and steps <ref>, <ref> and <ref> form one leapfrog step. Thus, the discretisation error is locally $O(h^3)$ [Neal, 2012]. Establishing the trajectory discretisation error is complicated, because it requires solving a partial differential equation (PDE). However, all found PDE solutions indicate that the error is $O(h^{2})$ [Leimkuhler and Matthews, 2015, Sivak et al., 2014], the same as HMC.

This scheme generalizes several popular samplers.
To obtain HMC, we may repeat BAB many times, and use O as momentum refreshment. Alternatively,
setting $\gamma = 0$ implies $a=1$, so the O steps disappear, and the integrator becomes one leapfrog step. Adding a separate momentum refreshment then also yields HMC. To obtain MALA[Leimkuhler and Matthews, 2015 recommend the BAOAB ordering for more accurate sampling in MALA.] or SGLD [Welling and Teh, 2011], we take $\gamma \to \infty$. This implies $a \to 0$, which means that step <ref> entirely ignores the incoming momentum $\vm_n$ and just samples a fresh one. Because of this, we do not need to calculate the next momentum $\vm_{n+1}$ at the end of a step, or represent the momentum at all. Thus, steps <ref> and <ref> fall out. We may also fold steps <ref> and <ref> into <ref> to obtain the more common single-step expression for SGLD.

Note that there are other modern methods that simulate Langevin dynamics while preserving reversibility. For instance, Matthews and Weare, 2018 propose a variant of ABOBA that corrects for the (known or estimated) covariance of the gradient noise, similar to SGHMC, to get the correct $O(h^2)$ statistics. Also, Zhang et al., 2020 develop an integrator that is asymmetric, but nevertheless realizable in the reverse direction, and thus has a nonzero acceptance probability.

Acceptance probability.
The key ingredient is the ratio of backward and
forward transition probabilities, which we derive in <ref>.
The resulting acceptance probability is

\begin{equation}
\log \mhaccept{\vtheta_{n+1}, \vm_{n+1}}{\vtheta_{n}, \vm_{n}} = -\frac{1}{T}\bra{U(\vtheta_{n+1}) - U(\vtheta_{n}) + K(\vm_{n+3/4}) - K(\vm_{n+1/4})}.
\label{eq:ggmc-accept}
\end{equation}

This is the same expression as in HMC or Langevin MC, were we to regard only the inner
Leapfrog step as the transition to accept or reject.
Unlike Hamiltonian dynamics ($\gamma = 0$), underdamped Langevin dynamics ($0 < \gamma < \infty$) are not in general time-reversible, even when simulated exactly [Ma et al., 2019].[They are still backwards-realizable, that is, the backwards trajectory is in the support of the transition kernel.] However, applying a
Metropolis-Hastings correction step forces the resulting Markov chain to be reversible.
It is interesting, then, that the acceptance probability depends only on the
accuracy of the reversible steps BAB, and not
at all on the friction $\gamma$, which controls how much the momentum is refreshed and thus how irreversible the Markov chain is.

Empirical properties. In Fig. <ref>, we compare GGMC and HMC for different learning rates. GGMC has similar potential and kinetic energy for all chains. HMC has higher acceptance probability than GGMC (likely because the underlying chain in GGMC is irreversible). The behaviour of HMC varies from chain to chain, suggesting that it is unstable for larger $\ell$.
The only learning rate for which the acceptance probability is close to 1 is $10^{-6}$, which is much lower than is used in typical practice (see Fig. <ref> in the appendix).

Overall, the qualitative behaviour of GGMC seems similar to SGHMC, but with the advantage that its accuracy can be monitored or corrected using the M-H acceptance probability. See <ref> for more experiments and details.

Comparison of different learning rates ($\ell\sim h$, see <ref>),
for GGMC and HMC applied to a ResNet-20 in the CIFAR-10 data set; for epochs numbered 1000–1010. Details and full training run in <ref>.

§ ACCEPTANCE PROBABILITIES FOR MULTIPLE AND STOCHASTIC STEPS

Already Duane et al., 1987 observed that the Hamiltonian used to specify the transition probability $g$ need not be the one used to calculate the acceptance probability. Only the latter will determine what distribution the sampler converges to. That is, the sampler (HMC, or GGMC) may use stochastic or incorrect gradients, so long as the MH acceptance probability is calculated with respect to the target $\Ptarg{\vtheta}^{1/T}e^{-K(\vm)/T}$.

Unfortunately, calculating the exact acceptance probability after every GGMC step defeats the point of stochastic gradients. Instead, we can defer the acceptance probability calculation by many steps. This allows us to, for example, follow 10 epochs of stochastic gradients, and then spend one epoch calculating the acceptance probability. The resulting scheme multiplies the acceptance probabilities of each step together [Zhang et al., 2020]. Starting at $\vs_n$, after $N$ steps of MCMC, we can accept the resulting sample $\vs_{n+N}$ with probability $\min\bra{1, \mhaccept{\vs_{n+N}}{\vs_n}}$, where

\begin{equation}
\mhaccept{\vs_{n+N}}{\vs_n} = \prod_{i=1}^N \mhaccept{\vs_{n+i}}{\vs_{n+i-1}} = \frac{\Pboltzmann{\vs_{n+N}}}{\Pboltzmann{\vs_{n}}} \prod_{i=1}^N \frac{\Ptrans{\vs_{n+i}}{\vs_{n+i-1}}}{\Ptrans{\vs_{n+i-1}}{\vs_{n+i}}}.
\label{eq:multi-step-mh}
\end{equation}

Due to telescopic cancellation, we need only evaluate the exact likelihood at the end of many steps. This can be combined with other algorithms <cit.> to reduce the rejection probability.
<cit.> propose this deferred MH scheme and show its correctness for any sampler with a homogeneous Markov chain.
We extend it to symmetrically inhomogeneous Markov chains in <ref>. Our extension gives zero acceptance probability to asymmetric chains like the popular cosine learning rate schedule [Zhang et al., 2019].

§ CONCLUSION

We have shown that SGHMC and similar popular samplers cannot be corrected using rejection, since their acceptance probabilities are zero.
We have then revisited GGMC, which solves this problem and generalizes SGLD and HMC.
Finally, we have also presented a way to compute the acceptance probabilities of GGMC across several steps, leading to a method for drawing exact posterior samples using only stochastic gradients.

It is time for the Bayesian deep learning community to let go of the Euler-Maruyama scheme. Instead, we should explore the breadth of integrators for Langevin dynamics which we have inherited from molecular dynamics and statistics. We have shown that correctness is achievable. Thus, we should  demand correctness from our samplers.

The authors would like to thank Andrew Foong, Samuel Power, and Austin Tripp for helpful discussions and comments on an earlier draft. AGA was supported by a UK Engineering and Physical Sciences Research Council studentship [1950008]. VF was supported by a PhD fellowship from the Swiss Data Science Center.

[Bussi and Parrinello, 2007]
Giovanni Bussi and Michele Parrinello.
Accurate sampling using Langevin dynamics.
Physical Review E, 75, 2007.

[Chen et al., 2014]
Tianqi Chen, Emily Fox, and Carlos Guestrin.
Stochastic gradient Hamiltonian Monte Carlo.
In Proceedings of the 31st International Conference on Machine
Learning (ICML 2014), volume 32 of Proceedings of Machine Learning
Research, pages 1683–1691, 2014.
URL <http://proceedings.mlr.press/v32/cheni14.html>.

[Coelho, 2017]
Luis Pedro Coelho.
Jug: Software for parallel reproducible computation in python.
Journal of Open Research Software., 5:0 30, 2017.

[Ding et al., 2014]
Nan Ding, Youhan Fang, Ryan Babbush, Changyou Chen, Robert D Skeel, and Hartmut
Bayesian sampling using stochastic gradient thermostats.
In Advances in Neural Information Processing Systems,
volume 27, pages 3203–3211, 2014.

[Duane et al., 1987]
Simon Duane, A.D. Kennedy, Brian J. Pendleton, and Duncan Roweth.
Hybrid Monte Carlo.
Physics Letters B, 1950 (2):0 216 – 222,

[Geyer, 1992]
Charles J Geyer.
Practical Markov chain Monte Carlo.
Statistical science, pages 473–483, 1992.

[Hastings, 1970]
W.K. Hastings.
Monte Carlo sampling methods using Markov chains and their
Biometrika, 570 (1):0 97–109, 1970.

[Horowitz, 1991]
Alan M. Horowitz.
A generalized guided Monte Carlo algorithm.
Physics Letters B, 2680 (2):0 247 – 252,
ISSN 0370-2693.

[Leimkuhler and Matthews, 2015]
Ben Leimkuhler and Charles Matthews.
Numerical Methods for Stochastic Molecular Dynamics,
chapter 7, pages 261–328.
Springer International Publishing, Cham, 2015.
ISBN 978-3-319-16375-8.

[Leimkuhler et al., 2019]
Benedict Leimkuhler, Charles Matthews, and Tiffany Vlaar.
Partitioned integrators for thermodynamic parameterization of neural
URL <http://arxiv.org/abs/1908.11843v2>.

[Ma et al., 2019]
Yi-An Ma, Emily Fox, and Tianqi Chen.
Irreversible samplers from jump and continuous markov processes.
Statistics and Computing, 29:0 177–202, 2019.
URL <https://arxiv.org/abs/1608.05973>.

[Matthews and Weare, 2018]
Charles Matthews and Jonathan Weare.
Langevin Markov chain Monte Carlo with stochastic gradients.
URL <http://arxiv.org/abs/1805.08863v2>.

[Metropolis et al., 1953]
Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H
Teller, and Edward Teller.
Equation of state calculations by fast computing machines.
The Journal of Chemical Physics, 210 (6):0
1087–1092, 1953.

[Neal, 1993]
Radford M. Neal.
Probabilistic inference using Markov chain Monte Carlo methods.
Technical Report CRG-TR-93-1, 1993.
URL <https://www.cs.toronto.edu/ radford/review.abstract.html>.

[Neal, 1996]
Radford M. Neal.
Bayesian learning for neural networks, volume 118.
Springer, 1996.

[Neal, 2012]
Radford M. Neal.
MCMC using Hamiltonian dynamics.
URL <http://arxiv.org/abs/1206.1901v1>.

[Paszke et al., 2019]
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory
Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban
Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan
Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith
Pytorch: An imperative style, high-performance deep learning library.
In Advances in Neural Information Processing Systems 32
(NeurIPS 2019), pages 8024–8035. 2019.

[Roberts and Rosenthal, 1998]
Gareth O. Roberts and Jeffrey S. Rosenthal.
Optimal scaling of discrete approximations to Langevin diffusions.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 600 (1):0 255–268, 1998.

[Sivak et al., 2014]
David A. Sivak, John D. Chodera, and Gavin E. Crooks.
Time step rescaling recovers continuous-time dynamical properties for
discrete-time Langevin integration of nonequilibrium systems.
The Journal of Physical Chemistry B, 118, 2014.

[Sohl-Dickstein et al., 2014]
Jascha Sohl-Dickstein, Mayur Mudigonda, and Michael DeWeese.
Hamiltonian monte carlo without detailed balance.
In Proceedings of the 31st International Conference on Machine
Learning (ICML 2014), volume 32 of Proceedings of Machine Learning
Research, pages 719–726, 2014.
URL <http://proceedings.mlr.press/v32/sohl-dickstein14.html>.

[Welling and Teh, 2011]
Max Welling and Yee Whye Teh.
Bayesian learning via stochastic gradient Langevin dynamics.
In Proceedings of the 28th International Conference on Machine
Learning (ICML 2011), pages 681–688, 2011.

[Wenzel et al., 2020]
Florian Wenzel, Kevin Roth, Bastiaan S. Veeling, Jakub Świątkowski,
Linh Tran, Stephan Mandt, Jasper Snoek, Tim Salimans, Rodolphe Jenatton, and
Sebastian Nowozin.
How good is the Bayes posterior in deep neural networks really?,
URL <http://arxiv.org/abs/2002.02405v1>.

[Wu et al., 2020]
Hao Wu, Jonas Köhler, and Frank Noé.
Stochastic normalizing flows.
URL <http://arxiv.org/abs/2002.06707v3>.

[Zhang et al., 2019]
Ruqi Zhang, Chunyuan Li, Jianyi Zhang, Changyou Chen, and Andrew Gordon Wilson.
Cyclical stochastic gradient MCMC for Bayesian deep learning.
URL <http://arxiv.org/abs/1902.03932v2>.

[Zhang et al., 2020]
Ruqi Zhang, A. Feder Cooper, and Christopher De Sa.
AMAGOLD: Amortized Metropolis adjustment for efficient stochastic
gradient MCMC.
URL <http://arxiv.org/abs/2003.00193v1>.

§ ACCEPTANCE PROBABILITY CALCULATION FOR GGMC

The procedure used by Bussi and Parrinello, 2007 is to start with the joint
distribution of the noises $\vepsilon,\vepsilon'$ and transform it to the
distribution over $\vtheta_{n+1},\vm_{n+1} \vbar \vtheta_{n},\vm_{n}$ using the
density transformation formula with the determinant of the Jacobian.

However, in this case, there is somewhat more illustrative way. Consider

\begin{equation}
\vr = \sqrt{(1-a)T} \mM^{-1/2} \vepsilon,\quad\text{ and thus }\quad \vr\sim \Normal{{\boldsymbol 0}, (1-a)T\mM},
\end{equation}

and similarly for $\vr'$ and $\vepsilon'$. Note that $\vr = \vm_{n+1/4} - \sqrt{a}\vm_{n}$, and $\vr' = \vm_{n+1} - \sqrt{a}\vm_{n+3/4}$.
By substituting $\vr,\vr'$ into <ref>, we can write $\vtheta_{n+1},\vm_{n+1}$ as

\begin{align}
\vtheta_{n+1} &= \vtheta_n + h\mM^{-1} \bra{\sqrt{a} \vm_{n} - \frac{h}{2}\gradat{_{n}} + \vr}, \\
\vm_{n+1} &=
a \vm_{n} - \frac{\sqrt{a}h}{2}\bra{\gradat{_{n}} + \gradat{_{n+1}}} + \sqrt{a}\vr + \vr'.
\end{align}

The reciprocal of the determinant of the Jacobian is thus

\begin{equation}
\det{\begin{matrix}
{\dd\vtheta_{n+1}}/{\dd\vr} & {\dd\vtheta_{n+1}}/{\dd\vr'} \\
{\dd\vm_{n+1}}/{\dd\vr} & {\dd\vm_{n+1}}/{\dd\vr'}
\end{matrix}}^{-1}
= \det{\begin{matrix} h\vM^{-1} & 0 \\ \sqrt{a} & \mI \end{matrix}}^{{-1}} = h^{-1}\det{\vM},
\end{equation}

which is a constant, so it will be eliminated when dividing by the backward transition density.

Putting this all together, the forward transition density is

\begin{equation}
\begin{aligned}
\Ptrans{\vtheta_{n+1},\vm_{n+1}}{\vtheta_{n},\vm_{n}} = h^{-1}\det{\vM} \,
&\Normal{\vm_{n+1/4} - \sqrt{a}\vm_{n} \mvbar {\boldsymbol 0}, (1-a)T\mM}\\
&\Normal{\vm_{n+1} - \sqrt{a}\vm_{n+3/4} \mvbar {\boldsymbol 0}, (1-a)T\mM}.
\end{aligned}
\end{equation}

By substituting  $\vtheta_{n+1}\leftrightarrow \vtheta_{n}$,
$\vm_{n}\rightarrow -\vm_{n+1}$,
$\vm_{n+1/4}\rightarrow -\vm_{n+3/4}$,
$\vm_{n+3/4}\rightarrow -\vm_{n+1/4}$,
$\vm_{n+1}\rightarrow -\vm_{n}$,
we can obtain the backward transition probability

\begin{equation}
\begin{aligned}
\Ptrans{\vtheta_{n},-\vm_{n}}{\vtheta_{n+1},-\vm_{n+1}} = h^{-1}\det{\vM} \,
&\Normal{-\vm_{n+3/4} + \sqrt{a}\vm_{n+1} \mvbar {\boldsymbol 0}, (1-a)T\mM} \\
&\Normal{-\vm_{n} + \sqrt{a}\vm_{n+1/4} \mvbar {\boldsymbol 0}, (1-a)T\mM}.
\end{aligned}
\end{equation}

Dividing them, the constant terms from the Gaussians and $h^{-1}\det{\vM}$ disappear, and we obtain

\begin{equation}
\begin{aligned}
\log  & \frac{\Ptrans{\vtheta_{n},-\vm_{n}}{\vtheta_{n+1},-\vm_{n+1}}}{\Ptrans{\vtheta_{n+1},\vm_{n+1}}{\vtheta_{n},\vm_{n}}}
= \\
& -\frac{1}{(1-a) T} \bigg(
K\bra{-\vm_{n+3/4} + \sqrt{a}\vm_{n+1}} + K\bra{-\vm_{n} + \sqrt{a}\vm_{n+1/4}} \\
& \hspace{2.5cm} - K\bra{ \vm_{n+1/4} - \sqrt{a}\vm_{n}  } - K\bra{ \vm_{n+1} - \sqrt{a}\vm_{n+3/4}  }
\bigg)
\end{aligned}
\end{equation}

where $K\bra{\vv} = \frac{1}{2}\vv^{\tp}\mM^{{-1}}\vv$ is the kinetic energy defined in <ref>. This simplifies greatly, to

\begin{equation}
\log  \frac{\Ptrans{\vtheta_{n},-\vm_{n}}{\vtheta_{n+1},-\vm_{n+1}}}{\Ptrans{\vtheta_{n+1},\vm_{n+1}}{\vtheta_{n},\vm_{n}}}
=-\frac{1}{T} \bra{ K\bra{\vm_{n}} - K\bra{\vm_{n+1}} + K\bra{\vm_{n+3/4}} - K\bra{\vm_{n+1/4}} }.
\label{eq:forward-backward-log-ratio}
\end{equation}

The log ratio of the Boltzmann densities (<ref>) is

\begin{equation}
\begin{aligned}
\log \frac{\Pboltzmann{\vtheta_{n+1}, \vm_{n+1}}}{\Pboltzmann{\vtheta_{n}, \vm_{n}}} &= -\frac{1}{T}\bra{U(\vtheta_{n+1}) - U(\vtheta_{n}) + K(\vm_{n+1}) - K(\vm_{n})}.
\label{eq:boltzmann-log-ratio}
\end{aligned}
\end{equation}

Now, to obtain the acceptance probability in <ref> we add eqs. (<ref>) and (<ref>), obtaining

\begin{equation}
\tag{\ref{eq:ggmc-accept}}
\log \mhaccept{\vtheta_{n+1}, \vm_{n+1}}{\vtheta_{n}, \vm_{n}} = -\frac{1}{T}\bra{U(\vtheta_{n+1}) - U(\vtheta_{n}) + K(\vm_{n+3/4}) - K(\vm_{n+1/4})}
\end{equation}

which is exactly the same as <ref>.

§ RELATIONSHIP BETWEEN GGMC AND GRADIENT DESCENT PARAMETERS

Following [Wenzel et al., 2020], we found it useful to express the
parameters of the GGMC sampler (momentum permanence $a$ and time step $h$) in terms of the parameters used most commonly in gradient descent deep learning (learning rate $\ell$ and momentum $\beta$). This way, we can use parameters that have been found to work well by the deep learning community.

Gradient descent first updates the momentum, and then updates the parameters; in a time-asymmetric way, unlike GGMC. Thus, to obtain the parameter equivalences, we must start from the <ref> step, combine it with the <ref> step of the next iteration, and end at the <ref> step. The equivalences then are as follows, where $N$ is the number of data points in the training set,

\begin{align}
\beta &= a = e^{-\gamma h} &  \gamma &= - \sqrt{\frac{N}{\ell}} \log \beta  \\
\ell &= N \, h^{2} & h &= \sqrt{\frac{\ell}{N}}. \\
\intertext{The equivalence between $\beta$ and $\gamma$ is different for the
symplectic Euler discretization. In that case it is}
\beta_{\text{SE}} &= 1 - \gamma h & \gamma &= (1 - \beta_{SE})\sqrt{\frac{N}{\ell}}.
\end{align}

§ MULTIPLE STEPS IN SYMMETRICALLY INHOMOGENEOUS MARKOV CHAINS

We can extend the proof for the correctness of deferred acceptance probabilities <cit.> to inhomogeneous Markov chains
by considering the following scheme. The scheme will force the inhomogeneity to be symmetric in time.

We augment the state of the Markov chain with $N$ random variables, $\{\zeta_i\}_{i=1}^N$, that will determine the randomness of the intermediate steps. Their target distributions have to be time-symmetric, i.e. $p_{\zeta_i}(\cdot) = p_{\zeta_{N-i+1}}(\cdot)$ for all $i$. Before an MCMC multi-step, we sample all the $\zeta_i$ from their target distribution. Then, we use them to provide all the randomness for the intermediate steps, so each step is just the result of a deterministic transformation:

\begin{equation}
\vs_{n+i} = f_i(\vs_{n+i-1}; \zeta_1, \dots, \zeta_i),\text{ for all }i=1,\dots,N.
\end{equation}

The functions $f_i$ have to be time-symmetric as well ($f_i = f_{N-i+1}$). Otherwise, the backward transition is not realizable using the same $\zeta_i$ as the forward one, so the acceptance probability is zero. If they are time-symmetric, the resulting acceptance probability is <ref>.

In this view, the cosine schedule [Zhang et al., 2019] has acceptance probability zero, even if a backward-realizable integrator is used. The cosine schedule is not a time-symmetric transformation: the schedule of learning rates is high at the beginning and low at the end of a cycle.

Wu et al., 2020 employ a similar construction to use MCMC samplers as normalizing flows.
<cit.> proposed the “partial gradient” construction: a variant of this for HMC and dataset minibatching in particular, motivated by the same operator-splitting view of integrators as the leapfrog.

§ EXPERIMENTAL DETAILS

We use the ResNet-20 architecture from [Wenzel et al., 2020] on the CIFAR-10 data set with data augmentation. We vary the learning rate as $\lbrace 10^{-6}, 10^{-4}, 10^{-3}, 10^{-2}, 10^{-1} \rbrace$, discarding the last one because of its unstable results. Training runs for 3000 epochs using stochastic gradients with minibatch size 125 (which is divisible by 50000, the number of data points). Every 10 epochs (a cycle) we evaluate the acceptance probability exactly and store a sample. <ref> shows one such cycle, and <ref> shows the whole training runs. We do not apply a Metropolis-Hastings correction, but merely examine it to determine the correctness of the sampler.

Comparison of different learning rates
for GGMC and HMC applied to a ResNet-20 in the CIFAR-10 data set; for 3000 epochs. Samples shown are taken every 10 epochs. Quantities as displayed are averaged over the full CIFAR-10 data set, but a random data augmentation transformation is applied to each example.

We implemented the experiments in PyTorch [Paszke et al., 2019] for GPU acceleration and gradients, using Jug [Coelho, 2017] for job scheduling.


# CT respiratory motion synthesis using joint supervised and adversarial
learning

Y-H Cao1, V Bourbonne1,2, F Lucia1,2, U Schick1,2, J Bert1,2, V Jaouen1,∗ and
D Visvikis1,∗ 1 LaTIM, UMR Inserm 1101, Université de Bretagne Occidentale,
IMT Atlantique, Brest, France  2 CHRU Brest University Hospital, Brest, France
∗ D. Visvikis and V. Jaouen share senior authorship  vincent.jaouen@imt-
atlantique.fr

###### Abstract

Objective: Four-dimensional computed tomography (4DCT) imaging consists in
reconstructing a CT acquisition into multiple phases to track internal organ
and tumor motion. It is commonly used in radiotherapy treatment planning to
establish planning target volumes. However, 4DCT increases protocol
complexity, may not align with patient breathing during treatment, and lead to
higher radiation delivery. Approach: In this study, we propose a deep
synthesis method to generate pseudo respiratory CT phases from static images
for motion-aware treatment planning. The model produces patient-specific
deformation vector fields (DVFs) by conditioning synthesis on external patient
surface-based estimation, mimicking respiratory monitoring devices. A key
methodological contribution is to encourage DVF realism through supervised DVF
training while using an adversarial term jointly not only on the warped image
but also on the magnitude of the DVF itself. This way, we avoid excessive
smoothness typically obtained through deep unsupervised learning, and
encourage correlations with the respiratory amplitude. Main results:
Performance is evaluated using real 4DCT acquisitions with smaller tumor
volumes than previously reported. Results demonstrate for the first time that
the generated pseudo-respiratory CT phases can capture organ and tumor motion
with similar accuracy to repeated 4DCT scans of the same patient. Mean inter-
scans tumor center-of-mass distances and Dice similarity coefficients were
$1.97$mm and $0.63$, respectively, for real 4DCT phases and $2.35$mm and
$0.71$ for synthetic phases, and compares favorably to a state-of-the-art
technique (RMSim). Significance: This study presents a deep image synthesis
method that addresses the limitations of conventional 4DCT by generating
pseudo-respiratory CT phases from static images. Although further studies are
needed to assess the dosimetric impact of the proposed method, this approach
has the potential to reduce radiation exposure in radiotherapy treatment
planning while maintaining accurate motion representation. Our training and
testing code can be found at https://github.com/cyiheng/Dynagan.

Keywords: Image synthesis, GAN, dynamic imaging, respiratory motion,
radiotherapy

##  1 Introduction

Four-dimensional computed tomography (4DCT) consists in reconstructing a 3D CT
acquisition into multiple time frames to follow internal motion [7]. It is
used routinely for radiotherapy (RT) treatment planning of lung cancer to
ensure dose coverage to the tumor while reducing margins for better sparing of
surrounding healthy tissues. A 4DCT acquisition is generally achieved using an
external motion synchronization device (such as a respiratory belt) to average
the CT signal of several phases -or gates- of the respiratory cycle, a process
known as respiratory gating [34]. According to RT terminology, the Internal
Target Volume (ITV) including all potential locations of the Gross Tumor
Volume (GTV) during breathing is obtained from the 4DCT scan. After the
addition of set-up margins, a Planning Target Volume (PTV) is used for the
treatment planning [45].

The benefits of using 4DCT acquisition over conventional CT for ITV
delineation have been reported in many studies [48, 17, 35]. However, it also
suffers from several drawbacks. Aside from practical issues such as cost and
additional protocol complexity, the longer acquisition times result in much
higher radiation exposure, up to six times a conventional 3D CT [49].

Figure 1:  Variability of 4DCT tumor tracking. Here, two patients underwent
multiple repeat 4DCT scans due to excess variability between corresponding
phases. Green, yellow, and red dots indicate the center of the tumor in scans
1, 2, and 3 respectively. Colored crosses show the corresponding tumor
position in the other scans. Images from CHRU Brest.

Another important limitation of 4DCT is that the actual respiratory pattern
during the RT session may differ due to the inherent irregular nature of
respiration. It is thus common to perform repeat 4DCT scans to assess this
variability and capture the full extent of tumor motion (with usually greater
confidence given to the second 4DCT scan due to a better familiarity of the
patient with the protocol [8, 47]). Repeat scans can even be performed more
than once if the tumor position between identical gates differs exceedingly
(Fig. 1), which is relatively common in the lower lobes [43, 50].
Nevertheless, despite these limitations, 4DCT is a key component of modern
lung RT due to improved accuracy of target delineation and better healthy
tissue sparing [46]. There is therefore a clear need to achieve similar
dosimetric objectives with more optimized protocols in terms of imaging-
induced dose and clinical workflow complexity.

In this work, we propose a method to generate patient-specific pseudo
respiratory motion synchronized chest CT images from static 3D CT scans and
external surface-based estimation using deep image synthesis. Contrary to DIR,
the objective is to synthesize motion from a single image (i.e. without
fixed/moving image pairs). To this end, we propose several key contributions:

  * •

At the architectural level, we generate patient-specific deformation vector
fields (DVFs) during breathing by conditioning a 3D U-Net encoder-decoder to
the respiratory amplitude, that can be obtained with any external respiratory
monitoring system. This is done through adaptive instance normalization of the
bottleneck of the encoder, a technique usually used for style transfer [24].
At test-time, the proposed model synthesizes any desired respiratory phase
from a static 3D CT image. This approach allows for the creation of dynamic,
time-varying images that can model motion of internal lesions and organs
during breathing.

  * •

Contrary to existing approaches derived from unsupervised deformable image
registration (DIR) that use reconstruction losses and explicit DVF smoothness
constraints during model training, we propose to supervise our model using
highly accurate lung DVFs and to jointly penalize the warped image and the DVF
magnitude through an additional adversarial objective. Doing so, we show
substantial performance improvements over more conventional losses, achieving
more faithful synthetic motion patterns.

  * •

We validate our approach extensively on real 4DCT images, not only by
measuring the global motion of the lungs using volumetric overlaps like
previous studies, but also by considering tumors and neighboring organs at
risk. In particular, we show that the agreement between synthesized and real
4DCT phases is similar to that obtained between two real repeat 4DCT
acquisitions of the same patient. Using a private dataset in addition to a
public one reported in previous studies, we also stress the importance of
considering smaller tumor populations for model validation.

  * •

For reproducibility and ease of use, we provide our training code and a
graphical application with pre-trained model weights that can produce dynamic
time series from any static lung 3D CT image. Our code is available at
https://github.com/cyiheng/Dynagan.

##  2 Related works

(a) Conventional DIR

(b) DL-based DIR

Figure 2: Comparison of typical appearances of deformation vector fields
obtained from state-of-the art (a) conventional lung DIR [51] and (b) deep
learning-based lung DIR [18]. Figure 3: Proposed convolutional DVF generating
architecture with scalar conditioning. The respiratory amplitude is injected
through an AdaIN mechanism at the bottleneck of the encoder. Given an input
image $X$ and a respiratory amplitude $\alpha$, the model learns a DVF
$\phi^{*}_{\alpha}$ from supervised labels $\phi_{\alpha}$. Then the input
image is warped with a spatial transformer to obtain the synthetic phase
$Y^{*}=X\circ\phi^{*}_{\alpha}$. An adversarial loss based on the
concatenation of the warped image and the magnitude of the DVF
$||\phi_{\alpha}||$ is used to improve the realism of the DVF without
explicitly constraining its smoothness.

Correlations between patient surface respiratory signal and internal organ
motion were reported in previous studies [44, 12, 10]. These correlations have
been used to build motion models for the estimation of internal organ motion
from surrogate data [36]. Different approaches to model motion have been
described in the literature, for example, using diaphragm displacement [4] for
tumor motion or skin surface data for the estimation of internal respiratory
motion [37]. Fayad _et al._ [11] generated synthetic 4DCT phases by
interpolating DVFs between end-of-inhale (EOI) and end-of-exhale (EOE) phases
obtained by DIR. Based on similar concepts, Huang _et al._ [25] predicted the
principal components of respiratory DVFs with neural networks using time-of-
flight surface information and a conventional lung DIR method [51]. Both these
works however require a 4DCT acquisition during deployment. Other works focus
on 3D CT to 4DCT image generation using unsupervised DL-based DIR between
training phases in order to apply the transformation to a static 3D image at
test-time. Chang _et al._ [5] proposed a method based on conditioning a deep
generative adversarial network (GAN) to a latent regressor using an
unsupervised DL-based DIR network to generate realistic 4D-XCAT phantoms for
simulation studies. Romaguera _et al._ [41] use similar concepts to generate
4D MRI volumes. As these methods use unsupervised DL-based DIR, they do not
require to know the actual spatial transformation to be learned between the
phases. Instead, similarly to conventional DIR methods that explicitly
optimize a parameterized transformation [51, 19, 28], unsupervised DL-based
DIR methods rely on maximizing (differentiable) similarity metrics while
constraining DVF smoothness. A major argument in favor of unsupervised
learning is that it is impossible to obtain real ground truth DVFs, and that
supervised training on synthetic transformations through random warping [9,
22] can lead to unsatisfactory results [14]. Unsupervised DL-based DIR,
however, has not yet reached the accuracy and robustness of conventional, non-
learning-based DIR [20, 21]. This is especially true beyond the cranial vault
(e.g. in the lungs or the pelvis) where larger anatomical deformations and
sharp intensity transitions may occur. As of today, no current state-of-the-
art DL-based DIR methods [18, 31] can achieve substantial accuracy in 4DCT
imaging, with a general tendency to produce overly smoothed DVFs, for
instance, leading to unrealistic modeling of motion at sliding interfaces
(Fig. 2). On the other hand, conventional methods such as [51] consistently
achieve sub-millimeter average target registration errors (TRE) across public
lung datasets even with large deformations [25, 28], which is better than
current learning-based techniques [18, 52, 6, 53]. For this reason, in this
paper, we propose to learn DVFs conditioned on the patient’s respiratory
amplitude in a supervised fashion leveraging state-of-the-art conventional
DIR.

In previous works, we have investigated the possibility of synthesizing
respiratory phases from static 3D images by conditioning an image-to-image GAN
architecture on 1D respiratory breathing traces [29, 2]. The overall motion of
the lung was reproduced, with a lowering of the diaphragm consistent with the
corresponding scalar amplitude injected into the model. The approach was
validated quantitatively using global metrics such as lung volume similarity
and global image quality metrics (mean absolute error - MAE, root mean square
error - RMSE and structural similarity - SSIM). However, the accuracy at
smaller scales was not satisfactory (e.g. little to no tumor motion), which
was not captured by the global evaluation metrics considered in the
experiments. Critically, tumor motion was not accurately reproduced,
especially in the lower lobes [3]. This renders the approach of little
practical value in RT treatment planning. More recently, similar global,
image-level validation strategies were considered in two unsupervised
learning-based approaches. Sang and Ruan performed deep interpolation between
EOE and EOI phases by generating spatio-temporally smooth DVF [42]. Lee _et
al._ [33] proposed a conditional seq2seq DVF learning architecture also using
explicit smoothness penalty, requiring only one 3D CT image and a respiratory
trace extracted from the position of the diaphragm. However, as previously
mentioned, unsupervised learning based on smoothness constraints on the DVF
may lead to poor realism, especially in the lungs [51], which is not captured
by the global evaluation metrics reported in these works. It is thus critical
for accurate tumor motion synthesis to further validate synthetic motion at a
more local level, supporting major changes in the methodology and the
validation presented in the current work.

##  3 Methods and Materials

Fig. 3 illustrates the proposed method’s overall workflow. During training,
the model learns a mapping between the EOE phase and the various phases of
real 4DCT images. Generation is conditioned on the input phase and an external
surface-based estimation $\alpha(t)$, which in clinical routine can be
extracted from different clinically used respiratory monitoring devices (e.g.
pneumatic belts, RPM, optical surface imaging). The generator maps to the
image-$\alpha$ input pair a 3-component DVF in a supervised fashion based
using a state-of-the-art conventional DIR method [51]. During training, the
output DVFs are compared to their ground truth following an $\ell_{1}$
minimization objective. A key contribution is to further improve
transformation realism through an adversarial term based not only on the
warped input, as is typically done in DL-based DIR, but also on the magnitude
of the generated DVF itself.

At test-time, given a 3D CT input image and corresponding $\alpha$ values, the
model reproduces synthetic respiratory phase images that mimic breathing at
the specified amplitudes.

###  3.1 Problem formulation and notations

Let ${X}(\bm{x})$ be a 3D CT image, where $\bm{x}\in\mathbb{R}^{3}$ is the
image domain. Let also $\mathcal{Y}(\bm{x})$ be a 4DCT image, i.e. a dynamic,
vector-valued series of $N_{g}$ gates

$\mathcal{Y}(\bm{x})=\left({Y}_{0}(\bm{x}),...,{Y}_{k}(\bm{x}),...,{Y}_{N_{g}}(\bm{x})\right),$
(1)

where $k\in\mathbb{N}$ is the gate number, with $k=0$ corresponding to the EOE
phase without loss of generality.

Our objective is to synthesize a vector-valued dynamic 4D image mimicking the
breathing dynamics of a 4DCT acquisition from $X(\bm{x})$, i.e. a set of $N$
3D images :

$\mathcal{Y}^{*}(\bm{x})=\left(Y_{0}^{*}(X,\alpha_{0}),...,Y_{k}^{*}(X,\alpha_{k}),...,Y_{N}^{*}(X,\alpha_{N})\right),$
(2)

where $\alpha_{k}$ is the relative variation of respiratory amplitude (e.g. in
mm) with respect to $X(\bm{x})$. Doing so, we make the simplification that
there is a one-to-one (bijective) correspondence between the true image
$Y_{k}$ and respiratory amplitude $\alpha_{k}$, i.e. for each $\alpha_{k}$
there can only be one 3D CT image associated to it.

To this end, during training, we learn from a set of real 4DCT acquisitions
$\left\\{\mathcal{Y}^{i}\right\\}$. For each training 4DCT image, we learn to
reproduce a set of DVFs:

$\left(\phi_{1},...,\phi_{k},...,\phi_{N_{g}}\right):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3},$
(3)

between the EOE phase $Y_{0}(\bm{x})$ and the remaining $N_{g}$ phases. To
learn this transformation, we use a supervised DL-DIR network in the form of a
U-Net shaped convolutional encoder-decoder with skip connections $G$, the last
layer of which is a three-channel layer, one for each DVF dimension. To
achieve amplitude-dependent image synthesis, $G$ is conditioned by a scalar
value $\alpha_{k}$ that measures the relative respiratory amplitude between
phases $Y_{0}$ and $Y_{k}$, i.e. $\phi^{*}_{\alpha_{k}}=G(Y_{0},\alpha_{k})$
are the generated DVFs during training.

During deployment, DVFs are generated for any desired $\alpha$ value from a 3D
CT image to obtain the requested image after warping
$Y^{*}(X,\alpha)=X\circ\phi^{*}_{\alpha}$, where ${X\circ\phi^{*}_{\alpha}}$
represents $X$ warped by $\phi^{*}_{\alpha}$.

###  3.2 Scalar conditioning layer

To condition image generation to the respiratory pattern of the patient, we
use an adaptive instance normalization (AdaIN) layer [24] to adapt the
features of $G$ to the respiratory amplitude $\alpha$ at the level of the
bottleneck of the encoding branch. The original AdaIN layer applies the style
of an image to a reference image in tasks such as style transfer by
normalizing feature statistics using a variable channel-wise mean and
variance. Here, the layer is modified to update the features based on $\alpha$
instead :

$\textrm{AdaIN}(\bm{x},\alpha)=\gamma(\alpha)\left(\frac{\bm{x}-\mu(\bm{x})}{\sigma(\bm{x})}\right)+\beta(\alpha),$
(4)

where $\gamma$ and $\beta$ depend on $\alpha$ after passing through a fully
connected layer and a one-dimensional convolutional layer.

###  3.3 DVF learning objective

The training loss of our DVF synthesis model is a dual supervised-unsupervised
compound loss :

$\mathcal{L}(G,D)=\lambda_{1}\mathcal{L}_{1}^{\phi_{\alpha}}(G)+\lambda_{2}\mathcal{L}_{\text{GAN}}^{Y,||\phi_{\alpha}||}(G,D),$
(5)

where $\lambda_{1}$ and $\lambda_{2}$ are parameters controlling the balance
between these two terms.

Departing from unsupervised methods that would typically constrain DVF
smoothness explicitly, we train $G$ in a supervised manner by adding a
$\ell_{1}$ reconstruction loss on the marginal components of the DVF:

$\mathcal{L}_{1}^{\phi_{\alpha}}=\sum_{xyz}{||\phi^{*}_{\alpha_{k}}-\phi_{k}||_{1}},$
(6)

where each component of the field is treated independently. $\phi_{k}$ is the
output of pTVreg [51], arguably one of the best lung DIR algorithm based on
isotropic total variation minimization, achieving consistently sub-millimeter
accuracy across several public lung datasets [25, 54].

To encourage the generation of realistic DVFs $\phi_{\alpha}^{*}$ without
explicitly modeling field smoothness, we propose to train the model under a
joint adversarial objective that concatenates both the warped image
$Y_{0}\circ\phi_{\alpha}^{*}$ and the magnitude of the DVF itself:

$\mathcal{L}_{\text{GAN}}^{Y,||\phi_{\alpha}||}=\mathbb{E}[\mathrm{log}D(X,Y_{k},||\phi_{k}||)]+\mathbb{E}[\mathrm{log}(1-D(X,Y^{*}_{k},||\phi^{*}_{\alpha_{k}}||)].$
(7)

The GAN loss evaluates the distribution alignment between the target image
$Y_{k}$ and the generated image $Y^{*}_{k}$, i.e. the EOE phase $X$ warped by
$\phi_{k}^{*}$. To further improve synthesis, the magnitude of the DVF is
included in the adversarial cost. The rationale for this is that the
information contained in the magnitude relates to the motion amplitude of the
patient’s respiratory pattern, which is associated to the scalar conditioning
value $\alpha$. Higher values of $\alpha$ should thus correlate with higher
magnitudes. We demonstrate in the experimental section that the addition of
this term is indeed key to achieving more faithful motion synthesis.

(a) Distribution of tumor volumes in the 4D-Lung (red) and CHRU (blue)
datasets.

(b) Metrics evaluated between EOE & EOI phases

Figure 4: Characteristics of the two 4DCT datasets used in our experiments.

##  4 Validation setup

###  4.1 Quantitative metrics

Different metrics were used for the validation of the proposed respiratory
motion synthesis model.

Reconstruction quality We first considered global reconstruction metrics MAE
and SSIM to evaluate the overall quality of the synthesized 4DCT phases
$Y^{*}_{k}$ compared to the ground truth phase $Y_{k}$.

Region-wise evaluation We then measured the variation in global lung volumes
by comparing both the synthetic and the true phase to the initial EOE phase
$X$. The lung volume variation (VV) was defined as:

$\textrm{VV}=\frac{|V_{B}-V_{A}|}{V_{A}}\times 100,$ (8)

where $V_{A}$ and $V_{B}$ are the lung volumes of the EOE phase and the target
phase respectively.

As already mentioned in the introduction, alongside global lung motion
assessment, it is also crucial for RT treatment planning to evaluate the model
accuracy at the tumor level, a problem somewhat overlooked by previous 4DCT
synthesis reports [42, 33]. To this end, we used the Dice Similarity
Coefficient (DSC) between the GTVs delineated on $Y^{*}_{k}$ and the true
phase $Y_{k}$. Considering the large variability of tumor volumes, we also
used the Euclidean distance between the Centers of Mass (CoM) of the tumor as
an additional volume-independent metric, hereafter referred to as the CoM
error.

Respiration affects not only the lungs but also peripheral organs at risk
(OARs), the localization of which is paramount for accurate RT treatment
planning. For this reason, we also evaluated volumetric DSC and CoM errors on
a per-organ basis in the heart, liver, stomach, spleen, kidneys, and pancreas.

DVF quality To assess the quality of the generated DVFs, we quantified the
total variation (TV) of the components of the DVF as in [51]. The rationale is
that well-behaved DVFs that are able to reproduce sharp transitions (e.g. at
sliding interfaces in thoracic imaging) show lower TV than exceedingly smooth
DVFs. Following common practice in DIR we also computed the negative
determinant of the Jacobian of the transformation $|J(\phi^{*})|$ to further
characterize the smoothness of $\phi^{*}$.

###  4.2 Datasets

We evaluated and validated the quality of our proposed dynamic image synthesis
model using two 4DCT datasets that show remarkable differences.

4D-Lung dataset We used a public dataset from The Cancer Imaging Archive
(TCIA), hereafter referred to as 4D-Lung [26] and reported in other 4DCT
synthesis studies [25, 42]. It is composed of $20$ images of locally advanced,
non-small cell lung cancer patients and includes GTV segmentation masks. Each
4DCT acquisition is composed of $10$ phases covering the full respiratory
cycle. The images were acquired on a 16-slice helical CT scanner (Brilliance
Big Bore, Philips Medical Systems) with a slice thickness of $3$mm and
$512\times 512$ axial resolution ($\sim 1$mm pixel size). Tumors are of rather
a large size with an average volume of $77\pm 96$ cm3. This translates into a
relatively high overlap in tumor volumes between the EOE and EOI phases,
assessed by a DSC of $0.72\pm 0.14$ and tumor CoM errors of $5.46\pm 2.66$ mm.

CHRU dataset Metrics reported on tumors such as found in the 4D-Lung dataset
hinder the fact that local accuracy may be insufficient to track smaller
volumes, due e.g. to the well-known dependency of the DSC on object size [40].
For this reason, we also considered a local clinical database acquired at the
radiation oncology department of Brest University hospital, France. This
dataset referred to as the CHRU dataset, contains patients with both primary
tumors and lung metastases. Due to quality concerns with certain 4DCT
acquisitions, not all data were included in the experiments. Images with
acquisition issues, such as missing thoracic portions or significant motion
artifacts, were excluded. Ultimately, only 43 patients were retained for the
experiments. Compared to the 4D-Lung dataset, the CHRU dataset shows much
smaller GTVs on average and greater volume variability, with mean volumes of
$6.4\pm 9.1$ cm3, which translates into much lower initial DSC $0.39\pm 0.31$
and greater distances between CoMs ($8.44\pm 6.04$ mm). Another notable
feature is inclusion of patients who underwent repeat 4DCT acquisitions, i.e.
30 patients had two 4DCT, five patients had three 4DCT, and one patient had
four 4DCT. This enables the assessment of inter-4DCT variability, which is a
key feature to evaluate the clinical interest of a respiratory motion
synthesis model. Each CT acquisition was reconstructed into eight respiratory
gates ($0\%$-to-$100\%$ inspiration followed by $75\%$-to-$25\%$ expiration by
steps of $25\%$). The images were acquired with a Siemens Sensation Open,
$2$mm slice thickness, and $512\times 512$ axial resolution ($\sim 1$mm pixel
size).

Fig. 4(a) shows the distribution of tumor volumes in the 4D-Lung and CHRU
datasets, with two typical examples. To further characterize these
differences, Fig. 4(b) shows the distribution of distances between the EOE and
EOI phases in each 4DCT acquisition, which relates to the amplitude of motion
correction required by the motion synthesis approach. Compared to 4D-Lung, the
CHRU dataset shows greater volumes differences between EOE and EOI phases,
smaller GTVs, and a greater diversity of tumor appearances.

(a) Reconstruction loss only ($\mathcal{L}_{1}^{\phi}$)

(b) w/ unsupervised GAN loss
($\mathcal{L}_{1}^{\phi}+\mathcal{L}_{\text{GAN}}^{Y}$)

(c) Proposed joint loss
($\mathcal{L}_{1}^{\phi}+\mathcal{L}_{\text{GAN}}^{Y,||\phi_{\alpha}||}$)

(d) Reference

Figure 5: Influence of the loss terms on the generated DVF. (a) Reconstruction
loss only: $\ell_{1}$ on the DVF (b) Addition of an adversarial term on the
warped image (c) Proposed compound loss: $\ell_{1}$ on the DVF + adversarial
term jointly on the warped image and the DVF magnitude. (d) Reference state-
of-the-art DVF obtained by [51]

###  4.3 Other experiments

4DCT variability study Due to the intrinsic variability of respiratory
patterns in patients acquired at different times, which relates to the
inherent limitations of 4DCT treatment planning, a useful metric to validate
the clinical relevance of the proposed respiratory motion synthesis method is
to ascertain that the distance between the synthetic phase $Y^{*}_{k}$ and
real phase $Y_{k}$ is not greater than the distance between repeat 4DCT
acquisitions $Y_{k}^{1}$ and $Y_{k}^{2}$ of the same patient. Repeat scans
were thus compared using DSC and CoM error to evaluate this variability.

Compound loss validation We quantified the improvement obtained by the
proposed compound loss based on supervised learning combined to an adversarial
loss on the magnitude of the DVF. To this end, we also trained alternative
models in which the total loss included combinations of reconstruction and
adversarial terms based on the warped image and/or the DVF.

Comparison to state of the art We compared our approach to a recent work,
RMSim [33], that also generates 4DCT sequences from static images and
respiratory traces. RMSim uses a modified seq2seq architecture based on
unsupervised DVF learning to enforce temporal consistency of the motion.
Contrary to our approach relying on external respiratory measurements, the
amplitude was obtained from image-derived internal information (maximum
displacement of the diaphragm in the superior-inferior axis), which is not
realistic and cannot be obtained in practice. The inference code was made
available by the authors of RMSim at the time of the completion of this work.
To ensure a fair comparison, we retrained our approach using similar image-
derived respiratory guidance, by measuring superior-inferior motion of the
diaphragm, since our method is adaptable to any form of scalar conditioning.
We did not retrain RMSim and used the network weights of the authors. We also
used the same training conditions for our method and learned a mapping from
EOI to EOE phase, rather than from EOE to EOI. The authors trained their model
on a private dataset using 100 patients with 10-phase 4DCTs using four A40
GPUs with 48GB memory each during 72 hours, which is considerably more than
our proposed model (see next section). The inference stage of RMSim required
35.2 GB of GPU memory. We compared the performances of these two models on
$20$ patients using the public 4D-Lung dataset.

###  4.4 Implementation and training

The model was trained in 3D using $42$ 4DCT images from $20$ patients from the
CHRU dataset. The testing set was composed by another $44$ 4DCT images from
$23$ additional patients of the CHRU dataset and $20$ 4DCT images from $20$
patients in the 4D-Lung dataset, which are also used for model comparison with
RMSim. The model was trained for $200$ epochs using ADAM optimization [32].
Loss hyperparameters were set to $\lambda_{1}=100$ and $\lambda_{2}=1$.

The initial 4DCT from the dataset were not aligned and had varying sizes in
the superior-inferior (SI) axis. We manually selected a reference 4DCT scan
with a suitable FOV along the SI axis and resampled all other 4DCT scans by
aligning the centers of origins. Images were then cropped around the lungs
using a robust DL-based automatic lung segmentation method [23]. We then
selected the largest bounding box (EOI phase), which we extended by $2$cm past
the lower point of the lungs to include abdominal organs of interest. Finally,
each image was resampled to a size of $128\times 128\times 128$ to allow the
neural network to be trained with relatively modest GPU capability. In all our
experiments, we used an NVIDIA GeForce GTX 1080Ti graphics card with 11GB of
memory and a batch size of 1.

Tumor masks available for each phase of the 4D-Lung dataset were used in
experiments requiring GTV delineations. For the CHRU dataset, GTVs were
obtained automatically in all phases using a nnU-Net model [27] trained with
the winning weights of the lung tumor task of the Medical Segmentation
Decathlon [1]. Heart segmentations were performed using PlatiPy, a free RT
segmentation software leveraging nnU-Net with additional image processing
heuristics [13]. Remaining organ segmentations (liver, stomach, kidneys, and
spleen) were obtained by training a nnU-Net model on the Multi-Modality
Abdominal Multi-Organ Segmentation Challenge (AMOS) 2022 training dataset
[30]. In all cases, final segmentation masks were reviewed by an expert and
manually corrected when necessary.

For each available repeat scan in the 4DCT variability study, we performed a
rigid registration between the EOE phases to avoid offsets in tumor and organ
positions. The CoM error was then computed between the two tumor
segmentations. The distance error in each axis was used to register the second
scan, correcting the initial CoM error. Subsequently, the transformation was
applied to the other phases and the segmentations of the repeat scan to ensure
consistency.

The scalar conditioning value $\alpha$ was computed from respiratory amplitude
measurements in a fully automated fashion. For each 4DCT phase of the two
datasets, we computed a respiratory amplitude (in mm) from the image following
the method described in [12], thus mimicking the knowledge that can be
obtained from a respiratory tracking belt such as the Anzai Respiratory Gating
System. Each phase was segmented to obtain the patient’s body [55]. Binary
difference images were obtained between the EOE and the other phases. The mean
distance in the anterior-posterior axis was defined as $\alpha$.

##  5 Results

###  5.1 Overall image and DVF analysis

To illustrate qualitatively the effect of the proposed compound loss used in
the training of our model, Fig. 5 compares generated DVFs obtained using other
variants of the total loss. When using only the reconstruction loss on the DVF
(Fig. 5(a)), the generated DVF was smoother and largely unidirectional along
the superior-inferior direction, with a low absolute magnitude corresponding
to a weaker lowering of the diaphragm. When an adversarial term was added only
considering the warped image $Y^{*}=X\circ\phi^{*}_{k}$, as can be done for
instance in unsupervised DL-DIR [15], the magnitude was closer to the
reference but nevertheless remained exceedingly smooth (Fig. 5(b)). More
faithful DVF appearance was achieved through the proposed dual supervised-
unsupervised loss taking also into account the magnitude of the DVF
$||\phi^{*}||$ in the joint adversarial cost (Fig. 5(c)). Sharper transitions
and orientations were achieved, in better agreement with the DVF obtained with
the reference registration method pTVreg (Fig. 5(d)).

To support these visual observations, Table 1 shows results obtained using
other combinations of loss terms. Training using only reconstruction losses in
the supervised setting $(\mathcal{L}_{1}^{\phi})$ or as in unsupervised DL-DIR
using the warped image $(\mathcal{L}_{1}^{Y})$ led to relatively high
reconstruction errors (mean MAE $\simeq$ 55 and mean SSIM $\simeq$ 0.57). When
the $\ell_{1}$ loss was combined to a GAN loss on the warped image $Y^{*}$
only (such as in [15] for DL-DIR), excessive DVF smoothness could be captured
by high TV and Jacobian values. On the other hand, best quantitative scores
with respect to all metrics considered were achieved when the adversarial term
jointly considered the magnitude of $\phi^{*}$, with much higher average SSIM
(0.74$\pm$0.05) and lower MAE (31.75$\pm$5.23), as well as smoothness measures
comparable to the reference DVF obtained through pTVreg.

Table 1: Compound loss analysis. | Warped image | DVF
---|---|---
| MAE | SSIM | TV | $|J|$
No synthesis | 54.53$\pm$8.98 | 0.53$\pm$0.06 | N/A | N/A
$\mathcal{L}_{1}^{Y}$ | 55.77$\pm$13.17 | 0.57$\pm$0.10 | 0.25$\pm$0.08 | 1.037$\pm$0.018
$\mathcal{L}_{1}^{\phi}$ | 56.29$\pm$13.79 | 0.57$\pm$0.10 | 0.26$\pm$0.08 | 1.037$\pm$0.016
$\mathcal{L}_{1}^{Y}+\mathcal{L}_{\text{GAN}}^{Y}$ | 43.00$\pm$10.64 | 0.57$\pm$0.07 | 0.41$\pm$0.08 | 1.028$\pm$0.010
$\mathcal{L}_{\text{GAN}}^{Y,||\phi_{\alpha}||}$ | 42.59$\pm$13.85 | 0.63$\pm$0.10 | 0.21$\pm$0.02 | 1.004$\pm$0.002
$\mathcal{L}_{1}^{\phi}+\mathcal{L}_{\text{GAN}}^{Y}$ | 37.39$\pm$12.12 | 0.68$\pm$0.10 | 0.12$\pm$0.04 | 1.008$\pm$0.005
$\mathcal{L}_{1}^{\phi}+\mathcal{L}_{\text{GAN}}^{Y,||\phi_{\alpha}||}$ | 31.75$\pm$5.23 | 0.74$\pm$0.05 | 0.12$\pm$0.04 | 1.006$\pm$0.004
Reference DVF | N/A | N/A | 0.13$\pm$0.03 | 1.007$\pm$0.006

Table 2 summarizes reconstruction results assessed by quantitative metrics MAE
and lung volume variation (VV). In the 4D-Lung dataset, mean VV was $6.86\%$
and mean MAE was $35.71$ before motion synthesis, suggesting relatively
moderate volume changes between the EOE and the other phases. After motion
synthesis, VV was reduced to $1.50\%$ on average. In the CHRU dataset, higher
initial lung VV of $11.66\%$ was reduced down to $3.05\%$ on average after
motion synthesis. MAE was also improved substantially. Training our model on
CHRU only was not detrimental to the performances achieved on the 4D-Lung
dataset, suggesting good generalization ability even when datasets have quite
different properties.

Table 2: Overall evaluation analysis. | 4D-Lung | CHRU
---|---|---
Initial VV (%) | 6.86 $\pm$ 5.32 | 11.66 $\pm$ 5.13
After synthesis VV (%) | 1.50 $\pm$ 1.19 | 3.05 $\pm$ 2.26
Initial MAE | 35.71 $\pm$ 15.95 | 56.47 $\pm$ 10.45
After synthesis MAE | 29.87 $\pm$ 10.74 | 42.26 $\pm$ 6.35

###  5.2 Region-wise evaluation

(a) EOE phase

(b) Ground truth

(c) Ours

Figure 6: Global lung motion on the coronal and sagittal plane. The red and green lines represent the initial position of the diaphragm in the left lung and the tumor respectively. The green arrow shows the location of the tumor. Table 3: Inter-phase similarity metrics at the tumor level before (no motion synthesis) and after synthesis. | 4D-Lung | CHRU
---|---|---
Initial DSC | 0.85 $\pm$ 0.13 | 0.43 $\pm$ 0.30
After synthesis DSC | 0.89 $\pm$ 0.09 | 0.68 $\pm$ 0.22
Initial CoM error (mm) | 2.83 $\pm$ 2.49 | 6.93 $\pm$ 4.97
After synthesis CoM error (mm) | 1.52 $\pm$ 1.29 | 2.74 $\pm$ 1.69

Fig. 6 shows a synthetic 4DCT phase from a patient of the CHRU dataset with a
relatively small tumor ($0.54$ cm3) visible both in the sagittal and coronal
planes. The tumor and the diaphragm position after motion synthesis matched
well the ground truth image in both planes (Fig. 6(b)). 111Animated visuals
are available in Supplementary Files.

Fig. 7 shows box-plots of volumetric overlaps between EOI and EOE phases for
the organs considered in our experiments in both datasets. For all organs, DSC
was substantially improved after motion synthesis. This indicates that our
model is not only capable of modeling motion inside the lungs but can also
deform peripheral organs with good accuracy, a higher overlap after synthesis
in all organs (DSC $\geq 0.65$).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 7: Comparison of organs overlap between images for the EOI phase in
4D-Lung dataset (A) and CHRU dataset (B). (a) Lung (b) Heart (c) Liver (d)
Stomach (e) Spleen (f) Left Kidney (g) Right Kidney (h) Pancreas (i) Tumor.

Table 3 shows motion synthesis results at the tumor level in both datasets. As
previously discussed, the 4D-Lung dataset showed lower initial CoM errors
between tumors with a mean of $2.83$mm and a mean tumor DSC of $0.85$ due to
the fact that tumors are of rather large size and motion relatively. In
contrast, the CHRU dataset had larger mean initial CoM errors ($6.93$mm) and
lower mean initial DSC ($0.43$). After synthesis, the mean CoM error was
reduced to $1.52$mm for a mean DSC of $0.89$ in the 4D-lung dataset (resp.
mean CoM error of $2.74$mm and mean DSC of $0.68$ in the CHRU dataset).
Despite differences in tumor characteristics, the model was thus able to
displace tumors in both datasets in a realistic and accurate fashion.

Figure 8: Tumor motion synthesis accuracy along the respiratory phases in the
4DCT-Lung dataset (assessed by CoM errors in mm). Our approach is compared to
RMSim, a recent work based on unsupervised learning of DVF fields with
temporal regularization [33]. Table 4: 4DCT variability study results

| CoM error (mm) | DSC
---|---|---
| Tumor | Tumor | Heart | Liver | Stomach | Spleen | Right Kidney | Left Kidney | Pancreas
4DCT1 vs 4DCT2 | 1.97$\pm$1.05 | 0.63$\pm$0.25 | 0.90$\pm$0.04 | 0.88$\pm$0.05 | 0.67$\pm$0.13 | 0.70$\pm$0.09 | 0.69$\pm$0.14 | 0.68$\pm$0.20 | 0.58$\pm$0.10
Ours1 vs 4DCT1 | 2.35$\pm$0.08 | 0.71$\pm$0.15 | 0.91$\pm$0.02 | 0.91$\pm$0.03 | 0.78$\pm$0.05 | 0.78$\pm$0.08 | 0.72$\pm$0.15 | 0.67$\pm$0.15 | 0.69$\pm$0.11
Ours2 vs 4DCT2 | 2.86$\pm$1.32 | 0.67$\pm$0.25 | 0.93$\pm$0.02 | 0.90$\pm$0.04 | 0.79$\pm$0.07 | 0.78$\pm$0.07 | 0.67$\pm$0.16 | 0.69$\pm$0.17 | 0.63$\pm$0.15

Table 4 shows quantitative results corresponding to the 4DCT variability
study, where we evaluated the agreement between two acquired 4DCT scans
compared to the agreement between real and synthetic scans. Within this
context, the EOI phase of the first 4DCT acquisition (4DCT1) was compared
against the EOI phase of the second scan (4DCT2). Mean inter-scans tumor CoM
errors and DSC were $1.97$mm and $0.63$ respectively, stressing the inherent
variability in the observed tumor position between real repeat 4DCT
acquisitions. Comparing real and synthetic phases led to similar mean CoM
errors of $2.35$mm and $2.86$mm and mean DSC of $0.71$ and $0.67$ for the
tumor. Similar results were achieved regarding volumetric overlaps between
OARs, with a variability comparable between images of acquired repeat and
synthetic images.

Fig. 8 compares tumor synthesis accuracy in the 4D-Lung dataset between our
approach and RMSim [33]. Initial CoM errors for the EOI phase were $5.46$ mm
on average, suggesting moderate tumor motion. Our model consistently
outperformed RMSim across all respiratory phases, despite a smaller training
dataset and lower computational requirements during both training and testing.

##  6 Discussion

Our main contribution lies in training the proposed model using supervised
learning on highly realistic DVF labels while using an adversarial term
jointly not only on the warped image (i.e. similar to unsupervised DL-DIR
training) but also critically on the magnitude of the field itself, a method
that combines advantages of both supervised and unsupervised learning. The
motivation for this strategy is to improve DVF realism and to indirectly
encourage correlations with the respiratory amplitude that is passed
simultaneously to the model, without constraining the smoothness of the DVF
through explicit regularization (e.g. as in RMSim [33]). Our results
demonstrate that this indeed substantially improves reconstruction quality
compared to conventional reconstruction losses or when using an adversarial
term on the warped image only, such as is typically done in unsupervised DL-
DIR [15]. This result, although being original in medical imaging, agrees with
earlier observations that adversarial losses are more effective for image-to-
image translation when reconstruction is not overly constrained [16]. This is
highly relevant to our application and likely explains why our prior efforts
using direct reconstruction losses did not produce satisfactory results at
small scales [2, 3].

In RT, accurate estimation of the motion of the tumors and the neighboring
OARs is critical for calculating the dose delivered to the patient and thus
maximizing treatment efficacy and minimizing the risk of radiation-induced
side effects. Using only the output of necessarily imperfect DVFs to evaluate
model performance is not satisfying. For this reason, evaluation was performed
using ground truth 4DCT acquisitions. The method was extensively validated
using two datasets that show notable differences in terms of GTV
distributions. We evaluated the robustness of our results at multiple scales
considering global image fidelity, global lung motion, tumor motion, and also
motion of neighboring OARs. Unlike previous studies that only considered
global reconstruction metrics [42, 33], covering more aspects of the model’s
performance across different metrics and datasets allowed for a more
comprehensive evaluation of the properties of the model. To remove dependency
on object size, we calculated in addition to volumetric overlaps the Euclidean
distance between the CoM of tumors as a volume-independent metric. Our results
show that it is critical to not only rely on single volume-based metrics in
order to report tumor tracking quality more faithfully [25]. The reported DSC
values were lower between repeat 4DCT phases compared to errors synthetic
images, indicating that our the position of the tumors were closer. The COM
error focuses primarily on the tumor’s overall position, while the DSC
coefficient considers both the shape and size of the tumor. In the case of
repeated 4DCT scans, the tumor may undergo minor positional changes without
significantly altering its overall shape. This results in a relatively low COM
error, suggesting that the tumor’s motion is similar between scans. However,
the DSC coefficient may be lower if the tumor undergoes minor changes in shape
or size between scans. Our model’s predictions tend to closely resemble the
shape of the tumor in the initial 3DCT scan, which is then used to predict
subsequent phases. This similarity in shape between the initial 3DCT and
synthetic phases facilitates better overlap between the segmentations of
generated images and real 4DCT images, despite the slight positional
discrepancies attributed to the COM error.

Despite leveraging the respiratory signal and image to enhance the accuracy of
motion prediction, determining the actual tumor position remains a challenge
due to the variation of the breathing patterns and other factors that may
influence the tumor motion. A key result of our work is the demonstration that
the error associated with motion synthesis is comparable with mismatches
typically obtained between repeat 4DCT acquisitions of the same patient, which
relates to the inherent limitations of 4DCT-based treatment planning. The
average tumor CoM errors achieved were less than $3$mm, which is in line with
a recent study of [47] recommending PTV margins due to respiratory motion
between $3$mm and $5$mm depending on tumor localization. Our results can thus
be considered satisfactory in regard to the intrinsic variability of 4DCT
treatment planning.

While we used state-of-the-art DIR to supervise our model, we acknowledge that
the inaccuracies in the DVFs may lead to errors in image synthesis. This is a
typical problem in machine learning when using noisy labels for supervision.
However, we showed that these errors remain within acceptable limits when
compared to the natural variability observed in repeat 4DCT scans. Our model
being agnostic to the DIR method chosen, the model can be improved by
retraining using better DVFs, if available, to further refine the motion
representation and enhance its accuracy and precision. Nevertheless, we
recognize the importance of assessing the robustness of our method to
variations in the chosen registration algorithm, as well as developing active
mechanism to combat label noise. This aspect will be thoroughly investigated
in future studies to ensure that our approach remains resilient to changes in
registration methods.

Regarding respiration conditioning, we proposed to constrain the model with an
amplitude-based scalar value extracted from patient-specific external 1D
respiratory signal traces. To this end, we added a modified AdaIN layer to the
bottleneck of the encoder based on this value. Doing so, we assumed a one-to-
one correspondence between the amplitude and the 4DCT phase image, which
represents a simplification of the problem. The use of AdaIN allows to
consider more elaborate (e.g. higher dimensional) conditioning if available.
For instance, more complex respiratory patterns could be better captured
through clinically available real-time entire external surface measurements by
optical devices [38, 39]. Including such information in the conditioning may
thus further improve the model performance, which we will investigate in
future research.

##  7 Conclusion

In this work, we have proposed a method to generate patient-specific pseudo
respiratory motion with the clinical objective of reducing the use of 4DCT
imaging in RT treatment planning. The model generates arbitrary breathing
phases from static images by conditioning motion synthesis to a scalar value
extracted from an external surface-based estimation. We train our generative
model explicitly on state-of-the-art lung DVFs while encouraging motion
realism under an adversarial objective that jointly constrains the warped
image and the DVF magnitude.

The method was extensively validated using two 4DCT lung tumor datasets
showing different characteristics through a variety of images, organs at risk,
and tumor level quality metrics. Results show for the first time that errors
due to motion synthesis are not higher than differences between real repeat
4DCT acquisitions of the same patient, which is an important result towards
4DCT-free treatment planning. Future methodological developments will focus on
more complex respiratory conditioning while further validation will evaluate
the dosimetric impact of the proposed approach.

## References

## References

  * [1] M. Antonelli et al.  The Medical Segmentation Decathlon.  Nature Communications, 13(1):4128, July 2022.
  * [2] Y-H. Cao et al.  Image and volume conditioning for respiratory motion synthesis using GANs.  In IEEE NSS MIC 2021, Oct 2021.
  * [3] Y-H. Cao et al.  Patient-specific 4DCT respiratory motion synthesis using tumor-aware GANs.  In IEEE NSS MIC 2022, Oct 2022.
  * [4] Laura I Cerviño, Alvin K Y Chao, Ajay Sandhu, and Steve B Jiang.  The diaphragm as an anatomic surrogate for lung tumor motion.  Physics in Medicine & Biology, 54(11):3529, may 2009.
  * [5] Y. Chang et al.  A generative adversarial network (GAN)-based technique for synthesizing realistic respiratory motion in the extended cardiac-torso (XCAT) phantoms.  Phys. Med. Biol., 66(11):115018, May 2021.
  * [6] W. Chi, Z. Xiang, and F. Guo.  Few-shot learning for deformable image registration in 4DCT images.  The British Journal of Radiology, 95(1129):20210819, January 2022\.
  * [7] J. D. Cox, J. Y. Chang, and R. Komaki.  Image-Guided Radiotherapy of Lung Cancer.  September 2007.
  * [8] S. De Oliveira Duarte et al.  Use of 4D-CT for radiotherapy planning and reality in France: Data from a national survey.  Cancer/Radiothérapie, 23(5):395–400, September 2019.
  * [9] K. A. J. Eppenhof and J. P. W. Pluim.  Pulmonary CT Registration Through Supervised Learning With Convolutional Neural Networks.  IEEE Trans. Med. Imag., 38(5):1097–1105, May 2019.
  * [10] H. Fayad et al.  Patient specific respiratory motion modeling using a 3D patient’s external surface.  Med. Phys., 39(6Part1):3386–3395, 2012.
  * [11] H. Fayad et al.  Generation of 4-Dimensional CT Images Based on 4-Dimensional PET–Derived Motion Fields.  J. Nucl. Med., 54(4):631–638, April 2013.
  * [12] H. Fayad, T. Pan, J. François Clement, and D. Visvikis.  Technical Note: Correlation of respiratory motion between external patient surface and internal anatomical landmarks.  Med. Phys., 38(6Part1):3157–3164, May 2011.
  * [13] R. Finnegan et al.  Feasibility of multi-atlas cardiac segmentation from thoracic planning CT in a probabilistic framework.  Phys. Med. Biol., 64(8):085006, April 2019.
  * [14] Y. Fu et al.  Deep learning in medical image registration: a review.  Phys. Med. Biol., 65(20):20TR01, October 2020.
  * [15] Y. Fu et al.  LungRegNet: an unsupervised deformable image registration method for 4D-CT lung.  Med. Phys., 47(4):1763–1774, April 2020.
  * [16] R. Groenendijk et al.  On the Benefit of Adversarial Training for Monocular Depth Estimation.  Computer Vision and Image Understanding, 190:102848, January 2020\.
  * [17] M. Guckenberger et al.  Is a Single Respiratory Correlated 4D-CT Study Sufficient for Evaluation of Breathing Motion?  International Journal of Radiation Oncology, Biology, Physics, 67(5):1352–1359, April 2007.
  * [18] L. Hansen and M. P. Heinrich.  GraphRegNet: Deep Graph Regularisation Networks on Sparse Keypoints for Dense Registration of 3D Lung CTs.  IEEE Trans. Med. Imag., 40(9):2246–2257, September 2021.
  * [19] H. P. Heinrich, M. Jenkinson, M. Brady, and J. A. Schnabel.  MRF-Based Deformable Registration and Ventilation Estimation of Lung CT.  IEEE Trans. Med. Imag., 32(7):1239–1248, July 2013.
  * [20] M. P. Heinrich and L. Hansen.  Voxelmorph++ Going beyond the cranial vault with keypoint supervision and multi-channel instance optimisation.  Technical Report arXiv:2203.00046, arXiv, February 2022.
  * [21] A. Hering et al.  Learn2Reg: comprehensive multi-task medical image registration challenge, dataset and evaluation in the era of deep learning.  IEEE Trans. Med. Imag., pages 1–1, 2022.
  * [22] M. Hoffmann et al.  SynthMorph: Learning Contrast-Invariant Registration Without Acquired Images.  IEEE Trans. Med. Imag., 41(3):543–558, March 2022.
  * [23] J. Hofmanninger et al.  Automatic lung segmentation in routine imaging is primarily a data diversity problem, not a methodology problem.  European Radiology Experimental, 4(1):50, December 2020.
  * [24] X. Huang and S. Belongie.  Arbitrary style transfer in real-time with adaptive instance normalization.  In ICCV, 2017.
  * [25] Y. Huang et al.  Deep learning-based synthetization of real-time in-treatment 4D images using surface motion and pretreatment images: A proof-of-concept study.  Med. Phys., 2022.
  * [26] G. D. Hugo et al.  A longitudinal four-dimensional computed tomography and cone beam computed tomography dataset for image-guided radiation therapy research in lung cancer.  Medical physics, 44(2):762–771, 2017.
  * [27] F. Isensee et al.  nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation.  Nature Methods, 18(2):203–211, February 2021.
  * [28] V. Jaouen, P-H. Conze, G. Dardenne, J. Bert, and D. Visvikis.  Regularized directional representations for medical image registration.  Technical Report arXiv:2111.15509, arXiv, November 2021.
  * [29] V. Jaouen and D. Visvikis.  4D respiratory motion synchronized image synthesis from static CT images using GANs.  In IEEE NSS MIC 2019, page 3, Oct 2019.
  * [30] Y. Ji et al.  AMOS: A Large-Scale Abdominal Multi-Organ Benchmark for Versatile Medical Image Segmentation.  Technical Report arXiv:2206.08023, arXiv, September 2022.
  * [31] Z. Jiang et al.  A multi-scale framework with unsupervised joint training of convolutional neural networks for pulmonary deformable image registration.  Phys. Med. Biol., 65(1):015011, January 2020.
  * [32] D. P. Kingma and J. Ba.  Adam: A Method for Stochastic Optimization.  arXiv:1412.6980 [cs], January 2017.
  * [33] D. Lee et al.  RMSim: controlled respiratory motion simulation on static patient scans.  Phys. Med. Biol., 2023.
  * [34] J. Liu et al.  Evaluation of the combined use of two different respiratory monitoring systems for 4D CT simulation and gated treatment.  Journal of Applied Clinical Medical Physics, 19(5):666–675, 2018\.
  * [35] A. V. Louie et al.  Inter-observer and intra-observer reliability for lung cancer target volume delineation in the 4D-CT era.  Radiotherapy and Oncology: Journal of the European Society for Therapeutic Radiology and Oncology, 95(2):166–171, May 2010.
  * [36] J. R. McClelland, D. J. Hawkes, T. Schaeffter, and A. P. King.  Respiratory motion models: A review.  Medical Image Analysis, 17(1):19–42, 2013.
  * [37] J R McClelland, S Hughes, M Modat, A Qureshi, S Ahmad, D B Landau, S Ourselin, and D J Hawkes.  Inter-fraction variations in respiratory motion models.  Physics in Medicine & Biology, 56(1):251, dec 2010.
  * [38] S. Nazir et al.  Technical Note: Kinect V2 surface filtering during gantry motion for radiotherapy applications.  Med. Phys., 45(4):1400–1407, April 2018.
  * [39] M. Ranjbar et al.  Validation of a CT-based motion model with in-situ fluoroscopy for lung surface deformation estimation.  Phys. Med. Biol., 66(4):045035, February 2021.
  * [40] A. Reinke et al.  Common Limitations of Image Processing Metrics: A Picture Story.  Technical Report arXiv:2104.05642, arXiv, July 2022.
  * [41] L. V. Romaguera et al.  Probabilistic 4D predictive model from in-room surrogates using conditional generative networks for image-guided radiotherapy.  Med. Image Anal., 74:102250, December 2021.
  * [42] Y. Sang and D. Ruan.  A conditional registration network for continuous 4D respiratory motion synthesis.  Med. Phys., 2023.
  * [43] S. Sarudis et al.  Systematic evaluation of lung tumor motion using four-dimensional computed tomography.  Acta Oncologica, 56(4):525–530, April 2017.
  * [44] A. Schweikard, H. Shiomi, and J. Adler.  Respiration tracking in radiosurgery.  Med. Phys., 31(10):2738–2741, 2004.
  * [45] G. Sgouros et al.  ICRU REPORT 96, Dosimetry-Guided Radiopharmaceutical Therapy.  Journal of the ICRU, 21(1):1–212, December 2021.
  * [46] J. Thariat et al.  Past, present, and future of radiotherapy for the benefit of patients.  Nature Reviews Clinical Oncology, 10(1):52–60, January 2013.
  * [47] P. Trémolières et al.  Lung stereotactic body radiation therapy: personalized PTV margins according to tumor location and number of four-dimensional CT scans.  Radiation Oncology, 17(1):5, January 2022.
  * [48] R. W. M. Underberg et al.  Four-dimensional CT scans for treatment planning in stereotactic radiotherapy for stage I lung cancer.  International Journal of Radiation Oncology, Biology, Physics, 60(4):1283–1290, November 2004.
  * [49] J. R. van Sörnsen de Koste et al.  Renal mobility during uncoached quiet respiration: An analysis of 4DCT scans.  International Journal of Radiation Oncology, Biology, Physics, 64(3):799–803, March 2006.
  * [50] I. Vergalasova and J. Cai.  A modern review of the uncertainties in volumetric imaging of respiratory-induced target motion in lung radiotherapy.  Med. Phys., 47(10):e988–e1008, October 2020.
  * [51] V. Vishnevskiy et al.  Isotropic Total Variation Regularization of Displacements in Parametric Image Registration.  IEEE Trans. Med. Imag., 36(2):385–395, February 2017.
  * [52] D. Wang et al.  PLOSL: Population learning followed by one shot learning pulmonary image registration using tissue volume preserving and vesselness constraints.  Med. Image Anal., 79:102434, July 2022.
  * [53] J. M Wolterink et al.  Implicit neural representations for deformable image registration.  In Medical Imaging with Deep Learning 2022, 2022.
  * [54] H. Xiao et al.  A dual-supervised deformation estimation model (DDEM) for constructing ultra-quality 4D-MRI based on a commercial low-quality 4D-MRI for liver cancer radiation therapy.  Med. Phys., 49(5):3159–3170, May 2022.
  * [55] X. Zhou et al.  Automated segmentations of skin, soft-tissue, and skeleton, from torso CT images.  In J. Michael Fitzpatrick and Milan Sonka, editors, Medical Imaging 2004: Image Processing, volume 5370, pages 1634 – 1639. International Society for Optics and Photonics, 2004.