{'fraction_non_alphanumeric': 0.034937443091761375, 'fraction_numerical': 0.01487201969446599, 'mean_word_length': 4.667717889908257, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gaussian Processes for real-time 3D motion and uncertainty estimation during MR-guided radiotherapy

Respiratory motion during radiotherapy causes uncertainty in the tumor's location, which is typically addressed by an increased radiation area and a decreased dose. As a result, the treatments' efficacy is reduced. The recently proposed hybrid MR-linac scanner holds the promise to efficiently deal with such respiratory motion through real-time adaptive MR-guided radiotherapy (MRgRT). For MRgRT, motion-fields should be estimated from MR-data and the radiotherapy plan should be adapted in real-time according to the estimated motion-fields. All of this should be performed with a total latency of maximally 200 ms, including data acquisition and reconstruction. A measure of confidence in such estimated motion-fields is highly desirable, for instance to ensure the patient's safety in case of unexpected and undesirable motion. In this work, we propose a framework based on Gaussian Processes to infer 3D motion-fields and uncertainty maps in real-time from only three readouts of MR-data. We demonstrated an inference frame rate up to 69 Hz including data acquisition and reconstruction, thereby exploiting the limited amount of required MR-data. Additionally, we designed a rejection criterion based on the motion-field uncertainty maps to demonstrate the framework's potential for quality assurance. The framework was validated in silico and in vivo on healthy volunteer data (n=5) acquired using an MR-linac, thereby taking into account different breathing patterns and controlled bulk motion. Results indicate end-point-errors with a 75th percentile below 1mm in silico, and a correct detection of erroneous motion estimates with the rejection criterion. Altogether, the results show the potential of the framework for application in real-time MR-guided radiotherapy with an MR-linac.


Introduction
Motion during abdominal radiotherapy decreases the efficacy of treatments due to an uncertain tumor location. This uncertainty can be reduced in several ways, one of sponding radiation beam adjustments. MR-guided radiotherapy introduces several major technical challenges, one of which is to reconstruct accurate 3D motion-fields in realtime, from a stream of MR-data. The required speed for this reconstruction is determined by the expected velocities; for slowly-moving tumors such as prostate tumors 1 Hz could be sufficient, but tumors subject to respiratory motion requires at least 5 Hz (Keall et al., 2006;Murphy et al., 2002).
Here, we focus on the latter category, which requires MR-data acquisition and motion reconstruction with a frame rate of at least 5 Hz. Motion-fields can be estimated from MR-images by means of image registration (image-based), or directly from k-space data (k-spacebased) (Huttinga et al., 2020). Reconstructing 3D motionfields of the abdomen and thorax at this rate is currently still challenging. To increase the achievable frame rate of 3D motion reconstructions, prior assumptions of images and/or motion-fields need to be included in the reconstruction. A frequently used strategy is to exploit an a priori built model with a two-step approach: 1) a calibration or training phase to build a patient-specific motion or image model prior to radiation; 2) a real-time phase during treatment, in which 3D information is reconstructed from a minimal amount of rapidly acquired MR-data by exploiting the model from the training phase. Examples of imagebased methods proposed for radiotherapy include MR-SIGMA by Feng et al. (2020), and approaches based on interleaved orthogonal cine images (Bjerre et al., 2013;Mickevicius and Paulson, 2017;Paganelli et al., 2018;Stemkens et al., 2016). An example of a k-space-based method for radiotherapy is the authors' real-time low-rank MR-MOTUS (Huttinga et al., 2020(Huttinga et al., , 2021b. MR-SIGMA estimated 3D MR-images at 3.3 Hz, Stemkens et al. (2016) achieved 3D motion-field reconstructions at about 2 Hz using image registration and cine-MRI, and MR-MOTUS achieved 3D motion-field reconstructions directly from k-space data at 6.7 Hz. Although considerably different in the modelling aspect, all the methods mentioned above employed a twostep approach.
It should be noted, however, that this two-step approach relies on the assumption that the motion in the training and real-time phases are similar. Although this is likely true in most cases, several practical scenarios such as bulk motion, or a change in breathing pattern could reduce the validity of this assumption. This could therefore result in erroneous motion estimates, which -if left undetectedcould eventually lead to harmful radiation to the patient. To warrant the patient's safety in such scenarios, methods for real-time MRgRT should therefore ideally not only estimate motion in real-time, but should also provide some measure of reliability. In a practical setting, this could be used for real-time quality assurance during radiotherapy, e.g. to halt and resume the radiation treatment according to the degree of confidence.
In this work, we present a probabilistic framework to simultaneously quantify 3D motion and provide a measure of reliability in real-time, which thereby addresses two critical needs for real-time MRgRT. This work is an extension of the preliminary work presented as conference abstracts in Sbrizzi et al. (2019) and Huttinga et al. (2021a). Our framework is based on the previously discussed two-step reconstruction approach. Firstly, in the training phase, a model for respiratory-resolved motion-fields is built that allows to represent 3D motion-fields with few coefficients. Secondly, in the inference phase, these representation coefficients are estimated from three mutually orthogonal readouts of MR-data in real-time, thereby exploiting the motion model built in the training phase.
The idea to extract motion information directly from few readouts was motivated by the success of the authors' MR-MOTUS method (Huttinga et al., 2021c), which reconstructs low-dimensional motion-field representation coefficients from few readouts of k-space data. The idea to use three mutually orthogonal readouts in order to do this is motivated as follows. A single readout of k-space data that crosses the k-space center effectively contains a projection of the excited FOV in the readout direction. This can be seen by transforming the readout to image space. It therefore mostly contains information of motion in the direction of the readout. Consequently, a set of three mutually orthogonal readouts contains information of motion in all directions. Based on the two observations above, we hypothesized that the low-dimensional motion-field representation coefficients can directly be inferred from the three mutually orthogonal readouts. That is, we assume the motion-field representation coefficients are a function of the three readouts.
Here, we propose to learn this underlying function via a probabilistic machine learning regression technique called Gaussian Processes (GP). A GP requires a calibration phase to tune its internal parameters based on a small training set, which takes ≈ 0.5 seconds in this work. The trained GP can then be used for real-time inference of the posterior distribution of the 3D motion-field representation coefficients, given the three mutually orthogonal readouts. This step exploits the availability of a closedform analytical expression for the posterior of a GP, which allows for sub-millisecond computations (≈ 0.1 milliseconds per dynamic in this work). Combining the resulting motion-field representation coefficients with the motion model, this eventually allowed for 3D motion-field reconstruction at frame rates as high as 69 Hz. Moreover, the posterior distribution as inferred by the GP captures not only the most likely motion-field estimate corresponding to the input data, but also the corresponding estimation uncertainty. The latter provides a measure of both the model-related and measurement-related uncertainties (respectively the epistemic and aleatoric uncertainties). As a consequence, a measure of motion model reliability is provided. We therefore hypothesize that the GP posterior uncertainty can be used for real-time quality assurance, i.e. to detect potentially erroneous motion estimates of the proposed framework.
Altogether, the framework could perform simultaneous real-time 3D motion-field estimation with real-time quality assurance. We assessed the accuracy of the presented framework based on several in silico and in vivo tests, and tested our hypothesis regarding the value of the GP posterior uncertainty for real-time quality assurance. For the latter, we designed a rejection criterion based on the uncertainty that flags dynamics with potential erroneous motion estimates in real-time. Both the motion estimation and rejection criterion were evaluated on data simulated using a digital XCAT phantom and in vivo data of 5 healthy volunteers acquired on an MR-linac, thereby considering different types of breathing and bulk motion.

A general introduction to Gaussian Processes
A Gaussian Process (GP) (Rasmussen, 2003) models a Gaussian probability distribution over functions, and can therefore be considered as an extension of the multi-variate Gaussian distribution, which models a Gaussian probability distribution over vectors. Gaussian Processes are frequently applied to regression problems by assuming the following model for noisy measurements y t ∈ Y ⊂ R N Y at samples x t ∈ X ⊂ C N X : where ∈ X , ∼ N (0, σ 2 n I), and y : X → Y is the underlying process. In case a function is drawn from a GP, and evaluated at a finite collection of N T samples, x T := [x 1 , . . . , x N T ] ∈ R N X ·N T , the vertically concatenated corresponding function values y T := [y 1 ; . . . ; y N T ] ∈ R N Y ·N T follow a multi-variate Gaussian distribution. A GP is completely characterized by a mean function m(x), and kernel covariance function k(x i , x j ), which in turn specify the mean vector and covariance matrix of the corresponding multi-variate Gaussian distribution: Here, T is defined as the training set T : , and m T and K T ,T denote the mean vector and covariance matrix, which are computed by evaluating respectively the mean and kernel function for all training samples in T .
The kernel function characterizes the properties of the underlying process (in this case y), such as smoothness or periodicity. Without loss of generality, a zero-mean GP is typically assumed with m ≡ 0, in which case the GP is completely characterized by the kernel function k(x i , x j ). The kernel function is typically a function of the distance between its two inputs, and is parameterized by hyperparameters θ that determine its form: Evidently, correct tuning of θ directly influences the GP's regression performance. Therefore, two steps are followed to perform GP-based regression. In the first step, the GP hyperparameters θ are estimated by maximum likelihood estimation (MLE) on the joint Gaussian likelihood (2). The result is a fully-determined GP kernel function, specifying the properties of the functions that best fit the training measurements. In the second step, a posterior distribution over the function values y Q at samples Q is computed by conditioning on the training data T . This yields the Gaussian posterior distribution (Rasmussen (2003), Eq. [A.6]): with In particular, the diagonal elements of the covariance matrix Σ Q define the GP prediction uncertainty. One of the main strengths of the GP framework is the availability of the closed-form analytical expressions in Eq. (4)-(5) that completely characterize the posterior distribution in Eq. (3). The computations involve the ma- where N T denotes the number of training samples and N Y the number outputs.

Proposed framework: integrating Gaussian Processes
in motion modeling The technical challenge for this application is to perform a real-time reconstruction of a complete 3D motion-field and corresponding reconstruction uncertainties from few readouts that can be rapidly acquired. This work considers Gaussian Processes for this purpose for two main reasons. Firstly, for the targeted application the computations in Eq. (4)-(5) can be performed in the order of milliseconds. Secondly, the probabilistic nature of GPs gives access to reconstruction uncertainties, which will prove valuable for quality assurance.
We identify 3 practical challenges when trying to apply GPs in our context. Firstly, GPs suffer from the curse of dimensionality, and are therefore challenging to apply to high-dimensional inputs (e.g. > 300 total samples on the three mutually orthogonal readouts); distances in high-dimensional input-space become uninformative (Bengio et al., 2006), while these distances -as evaluated by the kernel function -form the basis of the GP theory.
Secondly, GPs are mostly applied to regression problems with scalar functions, since modeling correlations in outputs is a challenging task and requires a non-trivial extension of the framework (Bonilla et al., 2007). To overcome these first two challenges, we propose a linear compression of both the input space and the output space, and learn a GP that maps from the latent input space X ⊂ R N X to the latent output space Y ⊂ R N Y . The schematic pipeline from few readouts s t ∈ S ⊂ C N S , to motion-fields d t ∈ D ⊂ R 3N , then becomes Training phase Inference phase Step 2: Gaussian Process regression Φ Φ Step 3: Expansion with , consisting of the compressed spoke representations x i and motion-fields representation coefficients y i . In the inference phase three steps are followed: 1) SN spokes are compressed to xt = U H st; 2) the posterior distribution P (yt|T ) is computed with the GP from xt; 3) a 3D motion-field dt and a corresponding uncertainty map are inferred by expanding the GP posterior distribution with the motion model Φ.
where U and Φ denote orthogonal bases of the latent input and output space, and U H the Hermitian transpose of U. This pipeline is also visually outlined in Figure 1. Assuming N Y ≥ 1 uncorrelated dimensions, each element [y t ] j (j = 1, . . . , N Y ) of y t can be modelled with a separate GP, each of which takes x t as input. Now that the general pipeline is outlined, the third practical challenge can be identified. Training of the GPs in the pipeline (6) requires training sets T j = {(x t , [y t ] j )} t=1,...,N T , for j = 1, . . . , N Y . On the one hand, the GPs are required to constantly evaluate distances between x t 's. To do this consistently, all x t should thus measure the same k-space locations but at a different time instance, which in turn requires the same k-space locations for all t. On the other hand, the targets y t will be derived from 3D images, and thus require k-space samples from the entire range of k-space coordinates. To meet both requirements, we employed an interleaved acquisition of mutually orthogonal self-navigation k-space spokes (SN spokes) and golden-mean 3D radial (GM3DR) Johnson (2017). The SN spokes are acquired along the three mutually orthogonal axes: feet-head (FH), anterior-posterior (AP), and leftright (LR). Before further processing the k-space spokes are transformed to image space through an FFT along the readout direction. Next, the vertically concatenated SN spokes at dynamic t, denoted as s t ∈ S are processed to GP inputs x t ∈ X , and GM3DR spokes at dynamic t to GP outputs y t ∈ Y. For the SN spokes, the processing includes the computation of a compressed representation with respect to an orthonormal basis U. For GM3DR spokes, the processing comprises a respiratory-resolved image reconstruction, followed by an image registration, the construction of an orthogonal motion model Φ via PCA, and finally the computation of a representation with respect to Φ. More details on these processing steps will be given in the subsequent sections. The construction of the GP training sets T j from the acquired SN and GM3DR spokes is also visualized in Figure 1.
The inference pipeline Eq. (6) requires three main steps: 1) a compression of readout spokes to the latent input space using U, 2) regression between the latent input and output space by a GP with hyperparameters θ, 3) an expansion of the latent output space to 3D motion-fields using a motion model Φ. These steps are also visualized at the bottom of Figure 1. All three sets of parameters (U, θ, Φ) mentioned above are obtained from respiratoryresolved training data, which was in turn acquired by sorting all GM3DR spokes with a 1D respiratory motion surrogate extracted from the closest feet-head spokeFeng et al. (2016); Huttinga et al. (2021b). Details including how the required parameters U, θ and Φ in the pipeline Eq. (6) are obtained will be discussed next, organized in sections. The first step in the pipeline, i.e. a compression of the acquired space, will be discussed in Section 2.2.1. The second step, i.e. the GP regression, will be discussed in Section 2.2.2. The third step, i.e. the expansion from GP outputs to motion-fields, will be discussed in Section 2.2.3. Finally, we also investigate the potential of the framework outlined above for quality assurance. To this extent, we designed and evaluated a rejection criterion based on the GP's reconstruction uncertainty. The rejection criterion is evaluated in parallel with the inference pipeline, rejecting unreliable motion estimates due to undesirable motion. Details on this are discussed in Section 2.3.

2.2.1.
Step 1: Input space compression with U We empirically observed a sensitivity of the reconstruction to cardiac motion, which resulted in high-frequency oscillations on top of the respiratory motion. Since this is undesirable, the k-space readouts were first transformed to image space with an FFT along the readout direction, and subsequently samples on the readouts (radial spokes) were thresholded based on their temporal frequency spectrum; any sample with a large contribution from the cardiac frequency range was removed. This thresholding typically resulted in an automatic selection of the liver, and excluded the heart and large arteries. After thresholding, respiratory-sorting was performed on the spokes, and PCA was performed along the respiratory dimension on the respiratory-sorted matrix of spokes S := [s 1 , . . . , s 20 ], such that S ≈ UV H . Here, and in the rest of this work, the number of respiratory bins was set to 20. The s i (i = 1, . . . , 20) were obtained by averaging over all SN spokes that ended up in the i-th bin. The number of columns in U and V for the approximation S ≈ UV H was restricted to max(4, N Y + 1), where N Y is the rank of the motion model as described above. This restricts the compression of the inputs to no less than 4 samples. This was empirically observed to improve the flexibility of the framework to correctly perform the regression. Finally, data compression at time t was performed by x t = U H s t .

2.2.2.
Step 2: GP regression between latent input and output space For the GPs we chose to use a Matérn kernel with ν = 3/2 and automatic relevance determination (ARD). The Matérn class of kernel functions is an extension of the class of radial basis functions, and restricts the functions that can be drawn from the GP to be ν times differentiable. With ν = 3/2, this thus results in a GP that models once-differentiable functions. In general, kernel functions return the correlation between its two inputs. In case of an ARD kernel, the relevance of each input dimension in this correlation computation is specified by the length-scales l 1 , . . . , l N X . For the Matérn-3/2 kernel, we have Here, θ := {σ n , l 1 , . . . , l N X } denote the kernel hyperparameters, and we defined L := diag 1/l 2 1 , . . . , 1/l 2 N X and the norm of a vector a with respect to a positive definite matrix B as a B := √ a T Ba. Theoretically, the standard deviation σ n of the target noise in Eq.
(1) could also be optimized, but here we manually fixed σ n such that the 95% confidence interval approximately contained the neighboring bin's training targets (see first column in Figure 3). With inputs x t ∈ R N X , this resulted in a total of N X optimizable GP hyperparameters. Finally, N Y pair-wise training sets were constructed as T j := {(U H s t , [y t ] j )} t=1,...,20 , and N Y sets of hyperparameters -one for each GP -were obtained through MLE on Eq. (2). The total time of the training was about 0.5 seconds. For the GP hyperparameter optimization we made use of the Gaussian Processes for Machine Learning (GPML) Matlab toolbox (Rasmussen and Nickisch, 2010). Once the training is performed, the resulting GP can infer a posterior distribution P (y t |x t ) of the motion-field representation coefficients y t , given the compressed representation of the spokes x t = U H s t . This inference process took about 0.1 milliseconds per dynamic.

2.2.3.
Step 3: Expansion from latent output space to motion-fields To expand the posterior distribution on y t to motionfields, a motion model is required. To derive the motion model Φ, first respiratory-resolved MR-images I := [I 1 , . . . , I 20 ] ∈ C N ×20 , were reconstructed with an L1-ESPIRiT reconstruction without temporal regularization with the BART toolbox (Uecker et al., 2015), where N denotes the number of voxels. Subsequently, respiratoryresolved motion-fields were obtained by image registration (Zachiu et al., 2015) between the end-exhale dynamic and all other dynamics, and stored in a matrix, D = [d 1 , . . . , d 20 ] ∈ R 3N ×20 . Similarly as in the twophase motion-field reconstruction approaches, PCA was performed along the respiratory dimension of D to obtain a linear motion model Φ. Typically, most information of D is concentrated in its first principal components (Huttinga et al., 2021b;Li et al., 2011;Stemkens et al., 2016;Zhang et al., 2007), allowing to truncate the number of columns of Φ from 20 to N Y 20, without sacrificing too much representational power. Thus, Φ ∈ R 3N ×N Y , and D is approximately low-rank, with rank N Y . In all experiments in this work, we normalized all columns in Φ to unit norm, and determined the rank via the L-curve of the explained variance of D as a function of the number of components. This resulted in 1 ≤ N Y ≤ 3 for all performed tests, and a motion model which approximates high-dimensional mo- A posterior distribution of the motion-field at time t can be derived by applying the motion model to the posterior distribution of the representation coefficients -P (y t |s t ) := N (µ t , Σ t ) -as output by the GP. Using basic properties of the multivariate Gaussian distribution, we can derive The spatio-temporal uncertainties P t in the motion-fields can be obtained as the diagonal of the covariance matrix in Eq. (8), i.e. P t := diag(ΦΣ t Φ T ). Note that this uncertainty map should be treated with caution. For the derivation above it is assumed that Φ is valid, and that there no uncertainty associated with it. Hence, the spatiotemporal uncertainties in the motion-fields d t can only be reasonably interpreted when this assumption is true. Unlike the uncertainty on d t , the uncertainty in y t has a useful meaning also when Φ is no longer valid. This forms the basis of the rejection criterion, as will be discussed next.

Online quality assurance via posterior uncertainty
We propose a rejection criterion based on the posterior uncertainty, which flags dynamics with uncertain predictions as follows. A new measurement can be considered unreliable when the corresponding prediction uncertainty is substantially higher than those evaluated on the training set. To this extent, we defined the rejection threshold τ such that where σ t is a measure of uncertainty. Given this definition of τ , any prediction with an uncertainty exceeding τ will be rejected. Effectively, this ensures a false positive rejection rate -i.e. the chance of incorrect rejections -of approximately 1-α. Note that τ in Eq. (9) is equivalent to the α-th percentile over all σ t evaluated on the training set data. In all experiments in this work we set α = 0.95. Each of the N Y GPs outputs a measure of uncertainty (i.e. the variance), which are combined to a single measure of uncertainty in order to evaluate the rejection criterion. To this end, we took a weighted linear combination, where the weights were determined by the singular values in the motion model. This leads to where s j is the singular value corresponding to the j-th singular vector of the motion model Φ.

Data simulation and acquisition
For the numerical simulations, data were generated using the digital XCAT phantom (Segars et al., 2010) with manually added MR-contrast. XCAT requires two inputs, an abdominal and chest waveform, and returns dynamic volumetric images with respiratory motion (cardiac motion was disabled) and corresponding ground-truth motionfields. To simulate realistic motion, cos 4 waveforms were used, hysteresis was simulated through a phase-delay between the two input waveforms, and the variations in the extreme respiratory positions were randomly generated as 1% and 2% of the waveform amplitudes for respectively end-exhale and end-inhale. All simulated motion-fields were post-processed with cid-X  to ensure realistic and invertible motion-fields. A complex reference volume in end-exhale was generated by adding a smoothly varying phase, and the resulting motion-fields were applied to this reference volume by cubic interpolation, resulting in a series of dynamics. Multi-channel images were simulated by multiplying each dynamic with 8 (static) coil sensitivity maps simulated with BART (Uecker et al., 2015). Finally, 800 readouts of k-space data were simulated for each dynamic with a NUFFT (Barnett et al., 2019) evaluated on the GM3DR+SN trajectory, and complex noise was added to the resulting k-space data to achieve an SNR of 60. In total four types of breathing modes were simulated for our experiments: normal, chestonly, abdominal-only, and amplitude drifts. Chest-only and abdominal-only breathing were simulated by setting respectively the abdominal and chest waveforms to zero, amplitude drifts were simulated by adding a linearly increasing shift to the normal breathing waveforms. For each breathing mode, approximately five breathing cycles were generated, resulting in a total of 100 dynamics. The total data per breathing mode corresponds to an acquisition of approximately 6 minutes. As described later in Section 3.2.1, half of the normal breathing phase data were used for training, and the other half for testing. Hence, the in silico training was performed on data that could in practice be acquired in 3 minutes. This is similar to the in vivo experiments, as will be discussed later in Section 3.3.
For the in vivo acquisitions, data were acquired on a 1.5T Elekta Unity MR-linac (Elekta AB, Stockholm, Sweden) from healthy volunteers during free-breathing. All experiments were approved by the institutional review board, carried out in accordance with the relevant guidelines and regulations, and written informed consent was obtained from all volunteers prior to the experiments. In all experiments, we used the standard clinical 8-element radiolucent array with anterior and posterior coils and a steady-state spoiled gradient echo sequence (SPGR) with TR = 4.8 ms, TE = 1.8 ms, FA = 20°, FOV = 30 cm x 30 cm x 30 cm, BW = 540 Hz, resolution = 3 mm x 3 mm x 3 mm. We considered two modes for the data acquisition. In the first mode, data were acquired with a golden-mean 3D radial (GM3DR) kooshball trajectory (Chan et al., 2009), interleaved every 31 spokes with three self-navigation spokes (SN-spokes); along feet-head (FH), anterior-posterior (AP), and left-right (LR). This first mode is referred to as the GM3D+SN mode. Since the proposed framework can infer a 3D motion-field from only a single set of SN-spokes, this interleaved acquisition allows for inference at 1000/(31 × 4.8) ≈ 6.7 Hz. The second mode served to test the feasibility of high-speed inference. In this mode only the SN-spokes were acquired, which allows for inference at 1000/(3 × 4.8) ≈ 69 Hz. The second mode is referred to as the SN-only mode.

In silico experiments
3.2.1. End-point-error analysis and model applicability tests To validate our framework, we performed an in silico end-point-error (EPE) analysis. These EPEs were computed per voxel, as the magnitude of the difference between the reconstructedd t and ground-truth motion-fields d t : where [·] i denotes the i-th voxel. We considered the scenario where we train only on data simulated during normal breathing, and performed inference on all types of breathing: normal, chest-only and abdominal-only breathing, and amplitude drifts. The normal breathing data were separated in two sets, a train and a test set. We analyzed the results in two ways. Firstly, we compared boxplots of the EPEs of a model trained for normal breathing, and applied to all breathing patterns. We hereby separated the EPEs statistics in dynamics rejected by our framework and the non-rejected dynamics, which allowed to see the effect of the proposed rejection criterion on the EPEs. Secondly, we analyzed the correlation between the GP posterior uncertainty and EPEs.

In vivo robustness test
To assess the robustness of the proposed framework we trained the model on three minutes of free-breathing and performed inference over all remaining GM3DR spokes. Data were simulated with the digital XCAT phantom, and acquired from 5 volunteers on the MR-linac, as described in Section 3.1. The volunteers were not instructed regarding their breathing. The preparation phase was performed on three minutes of data, as described in Section 2.2.1, Section 2.2.2 and Section 2.2.3. Inference was performed on the remainder of the data according to the pipeline in Eq. (6), and finally the rejection criterion was evaluated as described in Section 2.3.
The observable motion in FH spokes is typically used as a surrogate for respiratory motion (Feng et al., 2016). The first principal component of the motion model most significantly contributes to the final reconstructed motion-field. Hence, the motion in the FH spokes should be similar to the motion in the first principal component of the motion model, i.e. the output of the first GP. To this extent, we qualitatively validated the output of the first GP (posterior mean) by comparison with the projection images obtained from these FH k-space spokes. These projections were obtained by performing an FFT along the readout direction of the spokes. The comparison was performed on data acquired about two minutes after the training data.

Feasibility tests of high-speed inference at 69 Hz
In order to perform the input data dimension reduction, train the GP, and build a motion model, the presented framework requires both golden-mean 3D radial kooshball and SN spokes in the training phase. In the inference phase, however, the golden-mean radial spokes are not required, and the use of only the feet-head spokes would lead to a much higher inference frequency. With this experiment we investigate the feasibility of such high-speed inference.
To this extent, data were acquired in two phases. In the first phase, GM3DR+SN data were acquired in order to train the model. In the second phase, only the SN spokes were acquired at 69 Hz. Data were continuously acquired from four healthy volunteers (volunteers 2-5) over 9 minutes, with the same acquisition parameters as described in Section 3.1. The first 7 minutes of data were acquired with GM3DR+SN, and the last 2 minutes with SN-only. To mimic a realistic setting, only the first three minutes of the GM3DR+SN were used to perform all preparation steps in the training phase. The rejection criterion was calibrated on the last 30 seconds of the 3 minutes of GM3DR+SN training data. Finally, inference was performed on all available SN-only data, which covered around 2 minutes for most volunteers. This inference was performed on very high temporal resolution. The motion model was built to represent average breathing motion, and it is therefore likely that the data acquired at this high temporal resolution will represent motion states that slightly differ from the average breathing motion states. This could decrease the correlation between the training data and the inference data, which could in turn result in a high uncertainty. To make the rejection criterion less sensitive to this effect, we evaluated it at 5 Hz rather than at 69 Hz by only flagging dynamics for which the current dynamic and all 13 preceding dynamics (acquired in around 200 ms) exceeded the rejection threshold. Note that this temporal resolution for the rejection criterion of 5 Hz is still sufficiently fast for MRgRT (Keall et al., 2006;Murphy et al., 2002). Alternatively, the SN input spokes could be averaged in time, allowing for a trade-off between SNR of the input data and the temporal resolution of the outputs. Preliminary results on this are shown in Figure 8.
We also evaluated the effectiveness of the rejection criterion based on the GP posterior uncertainty by instructing volunteer 3 and 5 in the final minutes of the SN-only phase. Volunteer 3 was instructed to perform in sequence: 1) normal breathing, 2) a switch to chest-only motion, 3) bulk motion, and 4) normal breathing. Volunteer 5 was instructed to perform 1) irregular breathing, and 2) a bulk motion of the abdomen/thorax. During these undesirable movements, we compared the rejected dynamics with another independent method based on the center-ofmass (COM). The 3D COM coordinates were computed from the 3D radial spokes (Anderson III et al., 2013), and the comparison was performed with the left-right and feethead COM coordinates from a single channel of k-space data.
Finally, we performed an additional analysis to gain further insight in the uncertainty estimated by the GP, and thereby the dynamics rejected by the proposed method. To this aim, we analyzed the temporal behavior of the total uncertainty σ t (see Section 2.3) over the course of the whole acquisition, which was around 10 minutes for some volunteers. A gradual build up of this uncertainty would indicate a gradual decrease in the correlation between real-time measurements and the training data, and could therefore potentially give insights in organ drifts that could affect the inference during the SN-only phase.

In silico experiments
The left part of Figure 2 shows the results of the correlation analysis between GP posterior uncertainty and EPEs, indicating a positive correlation. Additionally, a statistical correlation analysis was performed between all breathing patterns other than the normal breathing pattern that was used for training. For respectively abdominal-only, chestonly and amplitude drifts, this resulted in correlation coefficients ρ = 0.99±0.00, ρ = 0.99±0.00 and ρ = 0.90±0.02, which all indicate a strong positive correlation between the GP posterior uncertainty and EPEs.
The right part of Figure 2 shows notched boxplots of EPEs before and after application of the rejection criterion. The results for normal breathing data, which was not exactly the same as the training data, show an interquartile range (IQR) of 0.36 mm -0.78 mm. Here, the IQR is defined as the data between the 25 th and 75 th percentiles. For abdominal, chest and drifts, the IQR was respectively 0.36 mm -0.68 mm, 0.29 mm -0.61 mm and 0.43 mm -0.88 mm. Overall, the maximum EPE after rejection was 1.31 mm. The notches indicate the 95% confidence interval around the medians, and show strong statistical evidence that the EPEs of rejected dynamics are higher than non-rejected dynamics.
Both result in Figure 2 empirically confirm our hypothesis that the GP posterior uncertainty can be used for realtime quality assurance, since the uncertainty is highly correlated with the error in the motion estimates. Especially the right part of Figure 2 highlights that the rejection criterion based on the GP posterior uncertainty rejects estimates with high EPEs and preserve low EPEs during undesirable breathing patterns. Figure 3 shows the results of the training phase (first column), an overview of the inference phase (second column), and a zoom in a region with increased estimation uncertainty. The number of ranks in the motion model per volunteer -and thus number of GPs trained per volunteer -was determined through the L-curve approach on the explained variances, and resulted for respectively the digital phantom and volunteers 1-5 in the ranks R = 2, 3, 3, 3, 1, 3, 2. The figure shows the output of the first GP, since this is generally the best interpretable component. In the first column the confidence interval should be observed; this was manually set in the training phase to account for potential errors in the training targets. The results in the overview in column 2 show realistic motion traces, and differences in breathing pattern can be observed between the different volunteers (different rows). In the zoom in column 3 it can be observed that the uncertainty increases in several scenarios: during a switch from normal to abdominal-only breathing (for the digital phantom, t ≈ 25 s), during deep inhales (volunteers 4 and 5) and during deep exhales (volunteers 1 and 2). The increased uncertainty is especially pronounced for dynamics with breathing amplitudes outside the range of the training data, i.e. larger than average exhales or inhales. This observation is explainable, since the GP posterior uncertainty increases with the distance to the training data, as evaluated by the kernel function. Animated Figure 1 in the supporting files shows the posterior mean and spatial estimation uncertainty maps for volunteer 1, as derived in (8). The animated figure shows the inference over the first 35 seconds in the second column in Figure 3, and visualizes every 4 th dynamic with a total 60 frames at 4 Hz. Similar to the results in Figure 3, the respiratory traces appear smooth overall, however small high frequency oscillations can be observed which could be related to cardiac motion or measurement imperfections (eddy currents). We refer the reader to Appendix A for an overview of all animated figures.

In vivo robustness test
The validation by comparison with the projection images obtained from FH k-space spokes is shown in Figure 4. Moreover, the figure shows the results of the application of the rejection criterion, as discussed in Section 2.3. Overall the posterior mean of the first GP coincides with the projection images. The framework mostly rejected endexhale dynamics, which were likely outside the range of the training data. Interestingly, for volunteer 2 the frequency of rejections is relatively large at the beginning, and decreases after about 400 seconds in the acquisition. For volunteers 4 and 5, it can be observed that the framework also rejected deep inhales, which were substantially different from the average breathing in the training data. The change in the frequency of rejections could be an indication of a breathing pattern change. The large number of rejections in exhale could be due to drifts of internal organs, resulting in a different end-exhale position. These hypotheses are further analyzed in Figure 5.

Feasibility tests of high-speed inference at 69 Hz
In this experiment, inference was performed at 69 Hz, using only three mutually orthogonal readouts as input. The results are shown in Figure 6, which compares the GP inference with projection images (background). The Animated Figures 2-5 in the supporting files show the posterior mean and spatial estimation uncertainty maps for volunteers 2-5, as derived in (8). These animated figures show the inference over the first 35 seconds of the data in Figure 6. The animation shows every 40 th dynamic with a total 60 frames, visualized at 4 Hz. An exception to this is volunteer 4, for whom every 80 th dynamic is visualized at 4 Hz because of a low breathing frequency. We refer the reader to Appendix A for an overview of all animated figures.
Overall the inference during normal breathing returns plausible motion traces. Volunteer 2 shows a regular breathing pattern, which changes to larger breathing amplitudes halfway the scan, and eventually shallow breathing at the very end of the scan. Volunteer 3 and 4 show a regular breathing pattern, and volunteer 5 shows a very irregular breathing pattern. The irregularities in the breathing are mostly accompanied by larger uncertainty, which is most pronounced for motion states that fall outside the range of motion in the training phase. The spatial uncertainty map for volunteer 3 shows a larger uncertainty than the other volunteers, although the visualized uncertainty range of 0 mm -1 mm should be noted here. The relatively larger uncertainty could be related to the relatively large breathing amplitudes of volunteer 3.
For volunteer 2, hardly any dynamics were rejected until about 640 seconds in the acquisition, after which the breathing pattern seems to change to very shallow breathing, resulting in many rejections. For volunteers 3 and 5, hardly any dynamics were rejected before the instructed motion changes, which started at respectively 535 seconds and 718 seconds in the acquisition. For volunteer 4, almost all end-exhale dynamics were rejected. Since this occurs after more than 10 minutes in the acquisition, this could be caused by organ drifts that changes the internal positions of the organs in such a way that it leads to rejections. This hypothesis is further analyzed in Figure 5. The results show that for most volunteers the posterior uncertainty gradually increases over the course of the acquisition, which indicates a gradual decrease of the correlation between the acquired data and the data in the training set. This could be explained by physiological drifts, which are also visually observable in Figure 3 and the projection images over the course of the whole scan (not shown).
To evaluate the rejection criterion, volunteer 3 and 5 were asked to perform a specific motion that forced a deviation from the data in the training phase and would render the motion model inapplicable. The results are shown in Figure 7. The top figure shows a comparison with the projection images obtained with an FFT along the readout of the FH k-space spokes, as well as the evaluation of the rejection criterion, as discussed in Section 2.3. The bottom part of the figure shows the center-of-mass (COM) coordinates in left-right (orange) and feet-head (purple), extracted from the 4 th channel of k-space data. The COM coordinates serve as an independent visualization of the volunteer's behavior during this experiment. The black arrows indicate the changes in motion, which are accompanied by a change in the pattern of the COM. Volunteer 3 (left), changed breathing pattern at arrow #1, performed bulk motion at arrow #2, and returned to a nor- Figure 3. Qualitative results of GP motion and uncertainty estimations in a digital phantom and five volunteers. The blue dots represent the training targets, the red line the posterior mean as estimated by the first GP -which corresponds to the first principal component of the motion model -from three readouts, and the gray shaded areas are the 95% confidence intervals obtained from the estimation uncertainty as output by the GP. The first and second columns show an overview of respectively the training and the inference phase, the third column shows a zoom-in on a region with a slightly different motion pattern or increased estimation uncertainty. The increased uncertainty mostly occurs at the highest motion amplitude levels. Note that the units on the vertical axes are arbitrary; the scaling dependents on the magnitudes of the principal components in the motion model. However, in all cases the higher values indicate exhales, while the lower values indicate inhales.  mal breathing pattern at arrow #3. Volunteer 5 shows irregular breathing with a slight drift towards inhale up to arrow #4, and bulk motion after arrow #4. The framework successfully detected abnormal motion, and rejected many dynamics. It should be noted that Volunteer 3 appears to return to a different position after the bulk motion, since the COM in left-right is slightly higher in the final part of the figure. Although this is not reflected by the framework's rejections, in practice the large number of consecutive rejections preceding this event should already be indicative that the model's predictions are too uncertain, and that the motion model should thus be updated. For volunteer 5, almost all dynamics were rejected during bulk motion. Only two sets of consecutive dynamics were not rejected during bulk motion, during which the COM coordinates indicate a similar position as in the preceding normal breathing phase.

Summary of innovations
In this work we presented a probabilistic framework which addresses two major technical hurdles towards real-time adaptive MRgRT simultaneously: real-time 3D motion-field estimation and uncertainty quantification. The framework was built on the idea that low-dimensional motion information can be extracted from few readouts of k-space data Huttinga et al. (2021c). This idea was exploited via a two-step reconstruction approach, in which first a motion model was built, and subsequently the model's coefficients were inferred from the data.
For the inference, a probabilistic machine learning regression technique based on Gaussian Processes was used. Due to its probabilistic nature, this technique not only estimated the most likely motion model coefficients, but also provided a measure of estimation uncertainty. The inferred model coefficients combined with the motion model yielded the motion-fields required for MRgRT. The inferred uncertainty was hypothesized to be useful for real-time quality assurance during radiotherapy. This hypothesis was empirically confirmed in simulations, in which it was shown to enable the detection of erroneous motion estimates, which -if left undetected -could in practice result in harmful radiation to organs-at-risk.
As opposed to other machine learning techniques such as neural networks, the training of the GPs could be done in about half a second, and required only a minimal amount of training samples (20 in this work). Moreover, the inference time was very short (≈ 0.1 milliseconds per dynamic) due to the availability of a rapidly computable closed-form analytical expression for the GP's posterior distribution. Altogether, this makes the GP a natural fit in a real-time MRgRT workflow where little latency is essential.
The complete framework was extensively validated in silico with simulations using a digital phantom, and in vivo with MR-linac data of five volunteers, thereby taking into account different breathing patterns and bulk motions.

Potential impact
The presented framework can have several applications in MR-guided radiotherapy. Firstly, the framework allows to infer 3D motion-fields at 69 Hz. The inference at this speed is more than sufficient to resolve abdominothoracic motion during MR-guided radiotherapy (Keall et al., 2006;Murphy et al., 2002). Although 5 Hz should be sufficient for this application, the feasibility of inference at 69 Hz could open up possibilities of applying radiotherapy to tumors subject to cardiac motion, such as central lung tumors. Another application could be cardiac radio-ablation (Cuculich et al., 2017). This is an emerging non-invasive treatment technique of cardiac arrhythmias with highly focused radiotherapy, which could benefit from high-speed tracking of myocardial landmarks.
Secondly, the proposed framework not only estimates motion-fields, but also provides a measure of estimation confidence. We have demonstrated an example of how this measure of confidence can be used for real-time quality assurance by designing a rejection criterion based on this estimation confidence. In practice, this could be useful to detect unexpected motions during radiotherapy, such as a change of breathing pattern or bulk motions, which could result in erroneous motion estimates due to an unsuitable motion model built for normal breathing. Without any measure of quality assurance, such estimates could result in harmful radiation to organs-at-risk. With the real-time quality assurance proposed in this work, the treatments could be (temporarily) halted to assure the patient's safety during such potentially erroneous motion estimates.
Finally, the availability of time-resolved 3D motionfields over the course of a radiotherapy treatment could be used for retrospective dose accumulation calculations (Kontaxis et al., 2020). Such calculations provide insights in the actual dose deposited to the target tumor and surrounding organs-at-risk during a treatment, which can be taken into account to improve subsequent treatment planning.

Related work
Real-time 3D inference as proposed in this work was shown before. For example, our own method MR-MOTUS (Huttinga et al., 2021c) achieved 3D motion-field reconstruction at 6.7 Hz, (Li et al., 2010) reconstructed 3D CT volumes and motion-fields at about 4 Hz, MR-SIGMA (Feng et al., 2020) estimated 3D MRI volumes at 3.3 Hz, and cine-based MRI methods such as proposed in Stemkens et al. (2016) achieved 3D motion-field reconstruction at about 2 Hz. Slightly different type of methods are based on surrogate signal models (Andreychenko et al., 2017;Low et al., 2005;Mcclelland et al., 2013Mcclelland et al., , 2017Tran et al., 2020), and also inferred 3D motion-fields at high temporal resolution from (several) 1D surrogate signals. Recently, also several deep learning (DL) based methods were proposed for real-time 3D inference. For example, Romaguera et al. (2021) inferred 3D motion-fields from a pre-treatment volume and real-time 2D images, and also Vol. 4 Vol. 5 Figure 6. High-speed in-vivo inference and quality assurance. The posterior mean of the first GP (green), and rejection criterion evaluations (red), are compared with projections on the FH-axis for all four volunteers with SN-only data available (2-5). Volunteers 3 and 5 were instructed to perform motion which would render the motion model invalid. Volunteer 3 performed in sequence: 1) normal breathing; 2) a switch to chest-only breathing; 3) bulk motion; 4) normal breathing, and the first abnormal event started around 535 seconds in the acquisition, as indicated by the yellow arrow. Volunteer 5 performed in sequence: 1) normal breathing; 2) irregular breathing; 3) bulk motion, and the first abnormal event started around 718 seconds in the acquisition, as indicated by the yellow arrow. The GP posterior means are scaled to visually overlap the pattern in the projection images. Figure 7. This figure shows the result of the high-speed inference for volunteers 3 (left) and 5 (right), as discussed in Section 3.3.2 and Section 4.2.2. The bottom graphs show a comparison with center-of-mass (COM) coordinates in leftright (orange) and feet-head (purple). The black arrows, labeled 1-4, indicate four events where the motion was changed, as can also be observed from the changing patterns in the COM-coordinates; 1) volunteer 3 switches breathing pattern, 2) volunteer 3 performs bulk motion, 3) volunteer 3 returns to normal breathing, 4) volunteer 5 performs bulk motion.
predicted the motion for future timepoints. Terpstra et al. (2021) proposed TEMPEST, a network that estimates 3D motion-fields directly from highly undersampled, aliased, 3D images with a frame rate up to 5 Hz. Frequently, non-DL-based real-time inference methods exploit a low-rank motion model similar to the one employed in this work, which is typically obtained by a retrospective compression of motion-fields using principal component analysis (Cai et al., 2015;King et al., 2012;Mishra et al., 2014;Stemkens et al., 2016;Zhang et al., 2007), or by decoupling the motion-fields into spatial components and temporal components based on surrogate signals (Low et al., 2005;Mcclelland et al., 2013Mcclelland et al., , 2017. Our previous works (Huttinga et al., 2021b,c) also employ a low-rank model, but estimate its components directly from k-space data by solving a large-scale non-linear optimization problem (Huttinga et al., 2020).
Most of the methods above only estimate motion-fields, without any measure of confidence or motion model applicability. An exception is the work by King et al. (2012), where a low-rank 3D motion model was fit to incoming 2D navigator MR-images and the motion model's applicability was constantly evaluated as the image similarity after registration. A notable difference is that King et al. (2012) used 2D cine navigator images, whereas we used three 1D spokes. The network proposed by Romaguera et al. (2021) also outputs estimation uncertainties, but these were not used for quality assurance. The proposed rejection criterion is therefore most similar to the model applicability test in King et al. (2012).
The employed dimension reduction techniques for both the input and output space were proposed to make GPs fit in a high-dimensional pipeline. Similar dimension reduction techniques on the output space were previously proposed in PCA-GPs (Higdon et al., 2008), with the aim to overcome the same challenges of combining GPs with high-dimensional output data. In this work, GPs were used to simultaneously perform regression over multiple scalar functions with multi-dimensional inputs, thereby assuming no correlations between the individual GP outputs. Alternatively, correlation in the outputs could be modeled with multi-task GPs (Bonilla et al., 2007). Although this could improve the performance in theory, our preliminary results in (Sbrizzi et al., 2019) with multi-task GPs showed little to no improvement.

Points of improvement and future work
Several aspects of this work could be improved for a clinical application. The proposed GP framework provides all the utilities required for real-time adaptive MRgRT. Further improvements on the method are modular and not strictly necessary. We expect the largest gains in performance to reside in improvements to the motion model, for example by better image reconstructions, or better image registrations. The number of rejected dynamics could possibly be reduced whenever a more expressive motion model would be used, e.g. a model built on 3D+t time-resolved cine MR-images, as proposed in King et al. (2012).
We have extensively validated our framework in silico with end-point-errors and four different breathing patterns. Moreover, the performance was assessed in vivo by comparisons with FH projection images. For a practical application, more in vivo validations are required with e.g. manually tracked targets on cine MR-images. This will be subject of a future work.
Several volunteers indicate a sensitivity of the framework to high-frequency oscillations in the data. These oscillations could be physiological, e.g. due to local cardiac motion or related blood flows. In the current framework, these high-frequency oscillations could not be coupled to local motion since this was not incorporated in the motion model. More complex motion models could possibly allow to identify the source of these oscillations. Alternatively, since inference could be performed with a speed of 69 Hz -which is about 10 times the required speed for MRgRTthe oscillations could simply be filtered out by performing a temporal averaging on the GP inputs. The averaging would allow for a trade-off between temporal resolution and smoothed outputs. Preliminary results in Figure 8 show that this is indeed feasible.
In this work, only correlations in the data were considered, but temporal correlations were not yet taken into account. Such an extension would require to model temporal correlations, but could prove valuable for near-future predictions. A possible downside of including temporal correlations is that it could impose strong restrictions on the temporal behavior of the motion (e.g. periodicity), which is why we have not considered it in this work.
In general, it should be noted that all improvements discussed above are modular, and would require little to no changes in the general pipeline of the proposed framework.

Conclusion
We have presented a probabilistic framework for simultaneous real-time 3D motion and uncertainty estimation. The complete framework, including the rejection criterion, allowed to preserve low EPEs (75 th percentiles ≤ 0.88 mm) during four different breathing patterns in simulations. Without the proposed rejection criterion these breathing patterns would have resulted in EPEs up to almost 6 mm, which -if left undetected -could lead to harmful radiation to organs-at-risk. The framework estimated in vivo motion that corresponds well with FH projections. Moreover, it flagged dynamics during which bulk motion and changes of breathing patterns were performed. This flagging strategy could be used to ensure the safety of the patient by (temporarily) halting the radiation, and continuing whenever confidence in the estimated motion is restored. Altogether, the framework tackles two major technical challenges for real-time adaptive MR-guided radiotherapy, real-time 3D MR-based motion estimation and uncertainty quantification, and it thereby paves the way to reach the ultimate potential of the MR-linac. Figure 8. The effect of temporally averaging the input SN-spokes. Several reconstructions of the first GP for volunteer 3 are shown for which a symmetric moving mean filter is applied with a varying width: 1, 27, 40, 60, 80, and 200 dynamics. This averaging increases the temporal latency, as indicated in the titles, since the data required by the second half of the window lies in the future. However, evidently, the filtering also results in smoother reconstructions. Hence, temporal filtering allows to make a trade-off between the temporal resolution of the reconstructions and the SNR. Temporal averaging with a width of 27 dynamics would result in sufficient temporal resolution for real-time adaptive MRgRT (Keall et al., 2006;Murphy et al., 2002).

Appendix A. Animated figures
All supporting animated figures corresponding to this manuscript are available at https://surfdrive.surf. nl/files/index.php/s/scLts9nJYXfbLMx.
Animated Figure 1 This GIF shows the posterior mean and spatial estimation uncertainty maps for volunteer 1, as derived in Eq. (8). The animated figure shows the inference over the first 35 seconds in the second column in Figure 3, and visualizes every 4 th dynamic with a total of 60 frames at 4 Hz.

Animated Figures 2-5
These GIFs show the posterior mean and spatial estimation uncertainty maps for volunteers 2-5, as derived in Eq. (8). These animated figures show the inference over the first 35 seconds of the data in Figure 6. The animation shows every 40 th dynamic with a total of 60 frames, visualized at 4 Hz. An exception to this is volunteer 4, for whom every 80 th dynamic is visualized at 4 Hz because of a low breathing frequency.