Real-time MR-based 3D motion monitoring using raw k-space data

Authors: marius  krusen
floris  ernst
Venue: MIDL_io_2024_Conference
Published: 2024-01-31 17:58:20

Abstract: Due to its great soft-tissue contrast and non-invasive nature, magnetic resonance imaging (MRI) is uniquely qualified for motion monitoring during radiotherapy.
However, real-time capabilities are limited by its long acquisition times, particularly in 3D, and require highly undersampling k-space resulting in lower image resolution and image artifacts.
In this paper, we propose a simple recurrent neural network (RNN) architecture to continually estimate target motion from single k-space spokes.
By directly using the incoming k-space data, additional image reconstruction steps are avoided and less data is required between estimations achieving a latency of only a few milliseconds.
The 4D XCAT phantom was used to generate realistic data of the abdomen affected by respiratory and cardiac motion and a simulated lesion inserted into the liver acted as the target.
We show that using a Kooshball trajectory to sample 3D k-space gives superior results compared to a stack-of-stars (SoS) trajectory.
The RNN quickly learns the motion pattern and can give new motion estimations at a frequency of more than 230 Hz, demonstrating the feasibility of drastically improving latency of MR-based motion monitoring systems.

Main paper: 
 introduction: 
 In recent years, the development of MR-linacs (Lagendijk et al., 2014; Mutic and Dempsey, 2014) has opened up new avenues for great advancements in the field of radiation oncology. An MR-linac combines a linear accelerator with magnetic resonance imaging (MRI) to allow online adaptation of the treatment plan before and even during the treatment with a radiation beam, giving rise to MR-guided radiotherapy. Compared to computed tomography (CT), the predominant imaging modality in radiotherapy, MRI offers superior soft-tissue contrast and does not expose the patient to ionizing radiation (Pereira et al., 2014) . It also does not require a highly invasive fiducial marker implantation unlike many X-ray imaging-based methods (Tang et al., 2007; Fledelius et al., 2014; Bertholet et al., 2017) .One of the critical challenges in realizing the full potential of MR-guided radiotherapy is the fast and accurate monitoring of intrafraction motion. Due to respiration, cardiac motion and other factors, the position of the target volume can change significantly during the treatment. For the liver, motion amplitudes of up to 30 mm have been reported (Park et al., 2012) . This can lead to underdosage of the target volume and overexposure of surrounding healthy tissue. Tracking the target in real time can mitigate this problem, reducing the risk of complications and improving treatment outcome (Al-Ward et al., 2018) .However, acquiring high resolution 3D MRI data typically takes several minutes, which is too slow for real-time motion monitoring.Since many MR-based motion monitoring methods rely on MR images, they have to work on lower resolution or fewer images, like using interleaved 2D slices (Bjerre et al., 2013; Tryggestad et al., 2013) , or reconstruct images from highly undersampled k-space data (Mickevicius and Paulson, 2017) to achieve a high enough speed. This can be combined with other MRI acceleration techniques like parallel imaging or compressed sensing to further reduce acquisition times (Tsao and Kozerke, 2012) . Deep learning-based reconstruction methods have also been proposed to further improve the speed and quality of the reconstructed images (Hammernik et al., 2017; Muckley et al., 2021) . Few approaches work directly on k-space data, including k-space-based navigators (Kober et al., 2012) and methods using the phase correlation method to register measured data to some reference k-space data (Vaillant et al., 2014) . This is often combined with a patient-specific motion model built during a less time-critical training phase before the treatment (Liu et al., 2021; Keijnemans et al., 2022 ). Huttinga et al. (2020) were able to reduce the latency of their MR-MOTUS framework from around 5 minutes for 3D motion to 170 ms by also splitting the reconstruction problem into an offline preparation and an online inference phase (Huttinga et al., 2022) .Our approach is to sample the k-space with a radial trajectory and directly use the incoming radial k-space lines (spokes) to continually estimate the target motion using a recurrent neural network (RNN). This way, no additional image reconstruction steps are performed and new motion estimations can be given after the acquisition of each new spoke in just milliseconds. In a previous work, we demonstrated this approach on a 2D rigid head motion estimation task achieving sub-millimeter accuracy (Krusen and Ernst, 2024) . Motion was simulated by rigidly transforming slices of 25 MRI head scans and the k-space was sampled using a golden-angle radial trajectory. Additionally, it was shown that removing the peripheral values of each spoke increases generalization across different patients.Here, we extend this approach to a more realistic 3D motion estimation scenario. The 4D XCAT phantom (Segars et al., 2010) is used to generate realistic abdominal MRI data affected by respiratory and cardiac motion. Since XCAT is a numerical phantom, there are no interpolation inaccuracies in the generated volumes and k-space sampling can also be done without interpolation. This sampling of the 3D k-space is done using a kooshball trajectory (Chan et al., 2009) , the 3D generalization of the golden-angle trajectory, and compared to a stack-of-stars (SoS) trajectory (Chandarana et al., 2011) , which uses goldenangle sampling along the k x -k y plane but Cartesian sampling along the k z dimension. A slightly adapted RNN is trained on the generated data to track the motion of a spherical tumor inserted into the liver. 
motion and mri simulation: 
 The motion and image data of the abdominal region was mainly generated using the XCAT phantom (Segars et al., 2010) . The 4D XCAT phantom is a numerical phantom supporting a number of anatomic models of male and female bodies with thousands of anatomical structures that simulates cardiac and respiratory motion at any spatial resolution and offers both deformable vector fields and image volumes as output. For these MRI motion datasets the resolution was set to 1 mm × 1 mm × 1 mm resulting in volumes of size 320 × 320 × 320 centered around the liver. Ten different anatomic models were used to create a dataset and for each model ten sequences of 1 000 time steps were generated. The length of one time step was set to 4 ms to simulate a steady-state spoiled gradient echo sequence with a repetition time (TR) of 4 ms.For each sequence a spherical tumor with a random diameter between 10 and 20 mm was inserted at a different location inside the liver. To further increase variance between sequences, different heart and respiration rates were used, and heart and respiration rate variability were simulated to vary the length of cardiac and respiration cycles during a sequence as well. An average resting heart rate of around 65 beats per minute and an average free-breathing respiration rate of around 13 breaths per minute, as well as the variability values were taken from Schumann et al. (2021) and Kobayashi (2021) respectively. Additionally, the initial phase of the cardiac and respiratory cycle at the beginning of a sequence was randomized and the two XCAT parameters controlling the maximum diaphragm motion and the maximum expansion of the chest along the anteroposterior axis were varied after every breath.-3.787 8.462 -11.343 MR contrast and a low-pass filter are applied to the volumes before transforming them to k-space and adding noise. Finally, the k-space is sampled resulting in the spoke data for the dataset.Simulating k-space spoke data from the image volumes generated with the XCAT phantom involved multiple post-processing steps as shown in Figure 1 . First, MR contrast is applied to the volumes. The pixel intensities of different organs and tissues for this step were extracted from T1-weighted MR images from the Duke Liver Dataset (Macdonald et al., 2020) and randomly varied between sequences. Additionally, a low-pass filter was applied to blur the pixelated edges between tissues. The resulting MR image volumes were transformed to k-space using a fast Fourier transform and complex-valued white Gaussian noise was added to achieve a signal-to-noise ratio of 20. Finally, the k-space was sampled following a kooshball trajectory for one dataset and following a stack-of-stars trajectory for a second dataset. Both are radial trajectories that would require interpolation. Since the XCAT phantom allows generating the image volumes at arbitrary rotations, each volume was already rotated during generation according to the desired k-space trajectory to make sure that the sampled spokes will align with the grid of the image volumes and no interpolation is necessary. The result is a sequence of 1 000 k-space spokes of 319 sampling points each.The tumor motion simulated by the XCAT phantom by default is a translation along the sagittal axis (y-axis) and longitudinal axis (z-axis) based on the respiratory phase and parameters like the maximum diaphragm motion and the maximum expansion of the chest.To have motion in all three dimensions, a custom sinusoidal motion trajectory was used that resembles the default motion along the y-axis and z-axis but also includes a smaller motion along the x-axis. On average, the motion amplitude along the x-axis is 30% and along the y-axis is 50% of the amplitude along the z-axis. Cardiac motion is simulated and affects the overall image volumes and k-space data, but does not contribute to the motion of the tumor. The resulting motion data is a sequence of translational motion along all three axes.In addition to the two datasets for both trajectories with a TR of 4 ms, another one with a TR of 40 ms is generated for each, for a total of four datasets. The higher TR allows longer simulated sequences spanning multiple cardiac and respiratory cycles without generating much longer sequences in this computationally expensive process.GRU FC + Dropout 20% FC Output Preprocessing Network 319 × 1 ∈ C 31 × 2 ∈ R 16 × 2 ∈ [-1, 1] 128 64 3Figure 2 : The preprocessing steps and network architecture. A complex-valued input spoke is split into real and imaginary data, shortened to the central 10 %, scaled down and its two halves are averaged during preprocessing. The gated recurrent unit (GRU) layer learns the temporal relation of the data before the fully connected (FC) layers with ReLU activation and dropout in between give the motion estimation.Before passing through the network, the input data goes through a few preprocessing steps. First, the complex-valued spoke data is split into its real and imaginary parts, since most machine learning methods only work on real-valued input data. Then, peripheral data of the spoke is removed leaving only the central 10 %. This improves generalization and leads to a more stable and fast training, as previously shown (Krusen and Ernst, 2024) . To further improve generalization, the average of the two symmetric halves of the spoke is used. Finally, the data is scaled by a constant factor that is sufficiently large, so that all current and future values are most likely in the range of [-1, 1] .During training two methods are used to increase the amount and variety of training data. The first method is subsequence learning. For each sequence of 1 000 time steps in the dataset there are many subsequences that can also be used for training. In particular, all 999 subsequences that only start later than the full sequence provide distinct and new data, since the start condition changes, so both the progression of the internal state of the GRU is different and the target motion is relative to a different initial position. Out of those, the 500 subsequences that are at least 500 time steps long are used as additional training data.The second method is applying different augmentations to the training data in each epoch. This gives more variety to the data during training, especially since the subsequence learning introduces many repeat inputs of the same data which can quickly cause overfitting. Three different augmentations are applied to the input data every time before passing it into the network. First, each spoke is scaled by a random factor between 0.9 and 1.1. Since this only alters the magnitude, the resulting spoke is still a valid input for the resulting motion target. Second, each individual value is multiplied by a factor sampled from a normal distribution with a mean of 1 and a standard deviation of 0.1. This can be seen as additional noise that changes every iteration. Lastly, there is a 5 % chance for a spoke to be swapped with the next one in a sequence. While this slightly invalidates the target values, the time between subsequent spokes is short enough that the changes to the target values would be miniscule. Furthermore, a noise augmentation is also applied to the target data by adding a value sampled from a normal distribution with a mean of 0 and a standard deviation of 0.1.Network training is done using the Adam optimizer (Kingma and Ba, 2014) and a learning rate of 0.001 that decays with a factor of 0.75 every 5 epochs for 25 epochs. All subsequences of same length are processed in one batch resulting in a batch size of 80. The mean square error between estimated and actual motion at each time step is used as the loss function. The training was performed on an NVIDIA A100 using the Pytorch framework and took about one hour.The trained networks are evaluated using a 5-fold cross validation. Each fold consists of the 10 full sequences of two different body models, resulting in 20 sequences per fold. The mean absolute error (MAE) between the estimated and actual motion at each time step is used as the evaluation metric. The average MAE of all sequences and folds is then used to compare the performance of the networks.Furthermore, the jitter of the output signal is calculated to determine how smooth the output is. This is an important quality for actually compensating the estimated motion, since a jittery output could exceed the mechanical capabilities of a subsequent multileaf collimator (MLC) for example. The jitter is calculated as the first derivative of the error between the estimated and actual motion. The mean absolute jitter then quantifies how much the output signal oscillates around the target. To reduce the jitter and smoothen the output, a simple moving average filter can be applied to the output signal. A filter length of 3 implies delaying the output signal by 1 time step and is used here to evaluate the effect on the jitter and the MAE.The network trained on the Kooshball dataset is compared to the network trained on the stack-of-stars dataset for both a TR of 4 ms and 40 ms. To put their scores into context, two simple additional estimators are used. The first one is a constant estimator that always outputs the mean motion of the target over the whole sequence to show the best possible error if no motion estimation is performed. The second one is a hypothetical perfect estimator with latency, whose output is the exact target function but delayed by 200 ms. 200 ms is the maximum suggested latency for real-time motion monitoring (Murphy et al., 2002) and comparable to existing approaches like MR-MOTUS (Huttinga et al., 2022) . This estimator shows the impact of just latency on the error. The reported error and uncertainty that real approaches have is not considered here and would further increase the total error. 
results: 
 The inference time of the trained networks is around 0.3 ms on an NVIDIA GeForce GTX 970, making the TR the main factor for the latency of the approach. Assuming a TR of 4 ms, the network can give new motion estimations at a frequency of more than 230 Hz. The results of the cross validation are shown in Table 1 . The RNN trained on the Kooshball dataset has a much lower MAE than the RNN trained on the stack-of-stars dataset for both TRs. For the more realistic TR of 4 ms, the Kooshball RNN has an MAE of 0.56 mm, 0.80 mm and 1.86 mm along the x-axis, y-axis and z-axis, respectively, slightly outperforming the perfect estimator with a delay. Appendix A shows how the MAE varies greatly between sequences and during a respiration cycle, explaining the large confidence intervals. The dataset with a TR of 40 ms causes both RNN methods to perform worse, while the two additional estimators expectedly perform similarly. In particular, the SoS RNN got even worse than the constant estimator anymore, suggesting that it does not learn much anymore. Appendix B shows that the subsequence learning and using half spokes during preprocessing improve the accuracy of the RNN. Applying a simple moving average filter to the output of the Kooshball RNN reduces the jitter to 0.030 while only slightly increasing the MAE to 0.57 mm, 0.96 mm and 1.87 mm, respectively. The SoS RNN shows a similar result. Smoothening the output for the higher TR of 40 ms would also reduce the jitter, but due to the longer delay, the MAE would also see an undesirable increase. The sequences in Figure 3 show example outputs of both RNN methods for both TRs overlaid on the target motion. The Kooshball RNN is able to follow the target motion accurately but often underestimates the peaks of the target motion. For the SoS RNN this behaviour can also be seen, albeit less accurately and with a noticeable oscillation with a period of around 20 time steps. This is exactly the number of slices in a stack in the dataset, showing one of the problems of the SoS trajectory for continuous motion estimation. 
discussion: 
 This work presents an learning-based approach that can rapidly estimate 3D motion from raw k-space data. The 4D XCAT phantom was used to generate abdominal MRI data affected by realistic respiratory and cardiac motion. A spherical tumor is inserted into the liver to provide a tracking target for the motion estimation. As a first step and due to the limits of the XCAT phantom, the tumor does not rotate or deform during motion unlike some of the other generated organs. However, these effects are typically small compared to the overall displacement of a tumor (Xu et al., 2014) and the phantom offers a reproducible way to generate volume data together with the corresponding motion data at a temporal resolution that is not feasible with real data.The presented RNN learns the motion pattern of the tumor to quickly give accurate motion estimations. While the accuracy is not as high as in the rigid 2D case (Krusen and Ernst, 2024) , the MAE is still small and allows for smaller planning target volumes and shorter treatment times. A hypothetical perfect estimator with a latency of 200 ms already has a similar error just due to latency and it does not yet take into account the actual tracking error that existing comparable approaches have. to the accuracy could be achieved by simulating more varied motion, possibly in combination with a more complex network architecture. Due to the periodic nature of respiratory motion the RNN could learn the general respiration pattern and not taking variations between respirations into account, as example outputs in Figure 3 could suggest. Training on more varied and less periodic data worked well for the 2D case (Krusen and Ernst, 2024) and could not just improve accuracy but also allow accurate motion estimation of other non-periodic intrafraction motion types, like slow drifting motion or patient motion.The Kooshball trajectory was shown to be superior to the stack-of-stars trajectory for this task. As the 3D generalization of the golden-angle trajectory, it provides a faster and more even coverage of 3D k-space avoiding artifacts like the oscillation seen in the SoS output and resulting in a much lower MAE.The main advantage of the approach is its very low latency. The effect that even a delay of just 200 ms has, which is considered enough for real-time motion monitoring, can be seen in the result of the hypothetical perfect estimator. With a estimation frequency of more than 230 Hz, the approach can apply a simple moving average filter to smoothen the output without an impact on the MAE. Additionally, predictions into the future, which are necessary in a motion compensation system to account for the system latency, are much more accurate with a short latency. 
conclusion: 
 In this work, we present an RNN approach to continuously estimate 3D motion from raw k-space data. By directly using the k-space data coming in from a 3D radial trajectory, the network can estimate motion with minimal latency at a frequency of more than 230 Hz. To train and validate the approach, the 4D XCAT phantom was used to generate realistic abdominal MRI data with an inserted tumor affected by respiratory and cardiac motion. The Kooshball trajectory was shown to be superior to the SoS trajectory for this task, since it more evenly covers the entire 3D k-space. The speed of the approach allows smoothening the output with a simple moving average filter and demonstrates the feasibility of a very lowlatency MR-based motion monitoring system. Future work will focus on further improving the accuracy of the approach. similar overall MAE suggests. However, the variation in MAE between sequences is quite high, explaining the large confidence intervals of all methods in Table 1 . The large number of outliers of the RNN-based method seen in sequences 2 and 7 stems from uncertainty at the beginning of the sequences. In the initial time steps, when the RNN did not receive much data yet, the output can sometimes be less or even more accurate than during the rest of the sequence. Another cause of variation in error is the respiration phase. As seen in Figure 5 , the mean absolute error varies during a respiration cycle. The value for the current respiration phase used by the XCAT phantom is scaled between 0 and 1, with 0 being full exhale and 0.4 being full inhale. The RNN has the largest errors around the full inhale phase when the motion amplitudes is largest which is typical for motion estimation methods. The delayed perfect estimator on the other hand has much larger variations during a respiration cycle and the smallest errors around the full inhale phase. This is due to the speed of motion being the slowest around the full inhale and full exhale phases, whereas the speed is much faster and the errors much larger inbetween these phases. Since the perfect estimator is delayed it is more sensible to fast motion. every spoke. The resulting MAE is slightly higher than the presented RNN, especially for a TR of 4, confirming the advantage of this preprocessing step. The second experiment evaluates the effect of subsequence learning by training without using subsequences. To keep results somewhat comparable, the training without subsequences runs for 1000 epochs instead of 25 and the learning rate decays every 200 epochs instead of every 5. While this still results in less total sequences processed during training, the results are already pretty clear with no real improvement happening during most of the training. The MAE for these experiments is much higher than that of the presented RNN. The dataset with a TR of 4, which also has much less variation in the target motion, is particularly affected by this and leads to overfitted results worse than the constant estimator.  
 
 Appendix:  


Decision: This is an interesting work for rigid motion estimation for a spherical tumor from successive spokes of radial MR acquisitions to be used in MR-guided radiotherapy. The network is trained using data simulated from an XCAT phantom, and experiments are performed on data simulated using both a kooshball and stack-of-stars trajectory. Unfortunately I do not see much value of the results of this work rather than the developed RNN model since the MRI simulations are so far from being realistic. There are no coils involved (very important to estimate motion from a specific area from center of kspace) and the simulated data is not realistic for MRI for the purpose of motion correction. If the method is evaluated on any real MRI data including multi-channel kspace data and shown to work in such realistic scenarios, I believe it will add value and can be ready for publication.

Reviews: review: paper_summary: paper_summary: This paper develops an RNN that predicts the rigid 3D motion of a spherical tumor from successive spokes of radial MR acquisitions to be used in MR-guided radiotherapy.  The network is trained using data simulated from an MR phantom, and experiments are performed on data simulated using both a kooshball and stack-of-stars trajectory. Experiments on simulated data show that the proposed method yields lower motion estimate errors than a constant motion estimate and achieves results slightly worse than a registration-based estimator with a prohibitive inference time. Experiments further show that the RNN trained on data simulated with the kooshball trajectory outperforms the RNN trained on data simulated with the stack of stars trajectory.
main_review: main_review: The paper provides an intuitive, well-explained method, and the authors have made substantive edits to the paper including extensions to 3D motion, additional ablation studies that confirm the importance of a few components of the method, and improved visualization which makes it easier to see patterns in the error of the method. I thus recommend to accept the paper, but will keep my review at "Weak Accept" due to the highly simulated nature of the data that is used to test the method (perfectly spherical tumor, artificially constant intensity across organs except for uniform noise).
strength_weakness: strength_weakness: - All experiments are on simulated data (perfectly spherical tumors that are only moved by respiratory motion, without rotations) and it is unclear how well the results would translate to more realistic data
- The method only outperforms one extremely simple baseline (constant motion prediction) and does not quite match the more sophisticated registration-based baseline
- The paper is missing ablation studies for several components of the method
- The paper is missing more sophisticated analysis of the results that would help the reader understand when it fails and succeeds
questions: questions: - Was every pixel within a given organ was assigned the same intensity value?
- How do we know this method would extend to a more realistic setting with a non-spherical tumor where rotations matter?
- How much of a difference do the averaging of the symmetric spoke halves and the subsequence learning strategy make?
- In what circumstances does the observed method fail relative to the registration-based method?
paper_summary: paper_summary: This paper extends a previous works from 2d rigid head motions to realistic 3d trajectories of abdominal magnetic resonance imaging subject to cardiac and respiratory motion. The author leverage RNN to design an algorithm able to track a tumor in a liver from data recorded with a phantom, and show reasonable accuracy for potential use in the clinic.
main_review: main_review: Very good, I do not have any further comments.                                                                                                                                                            .
strength_weakness: strength_weakness: I did not see any major weaknesses in this work. Maybe one could expand a bit on how the present work extend the previous work of the author, which is only mentioned in the introduction rapidly. Also the application of this method is only on a single, home made dataset, and it could be interesting to apply it to other, known dataset with reference, or benchmark results. I am not versed enough in this topic to propose/point to some, but if this is possible, it could be a nice addition to showcase the accuracy of this work. I would see this as a bonus only.
questions: questions: I have no questions.
paper_summary: paper_summary: This work proposes to use a recurrent neural network (RNN) to predict motion during MRI acquisition directly from k-space. The sequence of radial k-space samples are used as timepoints, and a motion prediction (in terms of difference with t=0) is output at every step. Predicting directly from undersampled k-space saves time allowing theoretical real-time motion correction during radiotherapy, for example. The paper builds upon previous work from the authors by exploring new sampling schemes on 3D data, instead of 2D.
main_review: main_review: Considering initial comments and the subsequent discussion, the authors have done a tremendous job of addressing all concerns raised (for example on the similarity with prior work or the results obtained from this work) by extending their proposed work and reworking the text to clear up any confusion.
strength_weakness: strength_weakness: The paper feels very incremental with little added value over previous work. Even though an additional dimension is added to the dataset, the network still predicts 2D motion. The results are somewhat disappointing, unfortunately.
questions: questions: Please put forward the distinctions with prior work and explain the seemingly mild results.

Comments: title: Thanks, I have updated my score.
comment: Thanks for the detailed and thorough response to my and other reviewers' comments. In particular, I appreciate the effort to motion to 3D, the additional ablation studies, and the much-improved visualization.

One comment about the additional results -- in Figure 4, it looks like a large number of examples are marked as outliers for sequences 2 and 7. Please add a comment in the text explaining why this may be the case.

Overall, my major concerns have been addressed, although the simulation environment is still somewhat artificial (constant pixel intensity except for uniform noise within each organ, and perfectly spherical tumors). I have thus improved my score to Weak Accept.
> Thank you for your quick response. We understand that unfinished results can and should not impact the review. However, the data generation and training has since been completed, and we have revised the paper to include the results of the new datasets. We could confirm our preliminary results and the new datasets not only solidify the 3D aspect of our method but also improve its accuracy. We have elaborated on this a bit more in our general response to all reviewers.

Thank you for addressing this concern. The revised method and results, in my opinion, shine a brighter light onto the proposed work and truly distance it from prior work.

> The method being much faster than it needs to be is only partly true. The frequency of around 230 Hz of the motion estimation would be unaffected by a bigger model, that would make use of the those 4 ms. However, the latency between motion occurring and being estimated by our method would increase. We have changed the text to call the TR "the main factor for the latency" instead. Still, this very low latency leaves room for a more complex and expressive model without latency-induced errors becoming a larger problem and is something we want to explore in the future.

Thank you for correcting my misunderstanding of this aspect of the proposed method and adapting the text to prevent other readers from making the same mistake. 

Considering the two discussion points above, I have no further comments and will be revising my initial rating.
Thank you for your quick response.
We understand that unfinished results can and should not impact the review.
However, the data generation and training has since been completed, and we have revised the paper to include the results of the new datasets.
We could confirm our preliminary results and the new datasets not only solidify the 3D aspect of our method but also improve its accuracy.
We have elaborated on this a bit more in our general response to all reviewers.

The method being much faster than it needs to be is only partly true.
The *frequency* of around 230 Hz of the motion estimation would be unaffected by a bigger model, that would make use of the those 4 ms.
However, the *latency* between motion occurring and being estimated by our method would increase.
We have changed the text to call the TR "the *main* factor for the latency" instead.
Still, this very low latency leaves room for a more complex and expressive model without latency-induced errors becoming a larger problem and is something we want to explore in the future.
Dear reviewers,

thank you very much for your detailed reviews and helpful feedback.
In response to the simplicity of only estimating motion in two axes, we have overwritten the XCAT tumor motion trajectories to include periodic motion in the x-axis based on the respiratory phase as well, generated new datasets with this motion and retrained the networks.
The paper and its results have been updated to reflect this change.

While the general sentiment of the paper remains the same, the additional motion axis not only solidifies the 3D aspect of our method but also improves its accuracy.
Small changes to the mean amplitude of motion caused by the change lead to slightly smaller errors in the y-direction and slightly larger errors in the z-direction for the two benchmark methods (from an MAE of 1.07 to 0.98 in y, and from 1.80 to 1.96 in z for the perfect estimator at TR 4).
Our proposed method has a similar change in the z-direction (going from an MAE of 1.74 to 1.86), but sees a large improvement for the error in the y-direction, which was previously larger than the error in the z-direction despite a smaller average motion amplitude.
Now, the MAE in y is 0.95 instead of 1.80, which also better matches the scale of the motion amplitude.
As a result, the total error is now comparable to the error caused by just the latency of the perfect estimator, which does not yet take the actual tracking error into account that real methods have.
This highlights the advantage of our approach even in terms of accuracy, since existing methods have latencies in the range of the simulated perfect estimator and will do those additional latency-induced errors.

Besides this, we have addressed your other comments in our individual responses and updated the paper accordingly.
Thank you for your reply.

> The default tumor motion generated by the XCAT phantom is limited to two axes. However, we have now modified the tumor motion to also include periodic motion in the x-axis based on the respiratory phase. Since the data generation is a time-consuming process, we have so far only been able to generate datasets with this motion pattern for the Kooshball trajectory but not the stack-of-stars trajectory. The preliminary results of training on these datasets are not just comparable, but actually seem to improve the accuracy of the approach, at least in the y-direction. The data generation will not be finished in time for the rebuttal, but we would like to include these results in the camera-ready version of the paper if they do not change the general sentiment of the paper.

While the inclusion of this dataset would have been a great addition to the paper, cementing the 3D aspect of the proposed method (as I presume x motion would now also have been predicted), I unfortunately cannot base my review on future results.

> In regard to the "somewhat disappointing" results, we have clarified the perfect estimator comparison to put our results better into context. Considering the additional error that existing methods have just due to delay highlights the speed advantage of our method and should present the results in a more positive manner.

Unfortunately, as far as my comprehension of the paper goes, the takeaway from the results as presented remain similar. The proposed method using Kooshball trajectories remains quite worse than the perfect estimator with maximum delay, with much higher variance (as highlighted in Figure 4, looking at the number of outliers). The speed of the proposed method seems to be its main appeal (c.f section 4: "The main advantage of the approach is its very low latency."), however section 3 mentions "making the TR the limiting factor for the latency of the approach". The proposed method therefore seems to be much faster than it needs to be (0.3ms vs 4ms), and would maybe benefit from a bigger model, more expressive model ?
Dear reviewer,

thank you very much for your helpful feedback and your comments on the writing style!

We agree that the paper is incremental, especially the network architecture.
However, we believe that the step from 2D to 3D is a significant change.
Additionally, the target motion changed from being an entirely rigid transformation of the whole field of view to the motion of a tumor surrounded by deforming organs with potentially different motion patterns.
The default tumor motion generated by the XCAT phantom is limited to two axes.
However, we have now modified the tumor motion to also include periodic motion in the x-axis based on the respiratory phase.
Since the data generation is a time-consuming process, we have so far only been able to generate datasets with this motion pattern for the Kooshball trajectory but not the stack-of-stars trajectory.
The preliminary results of training on these datasets are not just comparable, but actually seem to improve the accuracy of the approach, at least in the y-direction.
The data generation will not be finished in time for the rebuttal, but we would like to include these results in the camera-ready version of the paper if they do not change the general sentiment of the paper.

In regard to the "somewhat disappointing" results, we have clarified the perfect estimator comparison to put our results better into context.
Considering the additional error that existing methods have just due to delay highlights the speed advantage of our method and should present the results in a more positive manner.
If we assume that the total error from latency and tracking is not entirely additive, we can conservatively estimate that the total error of algorithms like MR-MOTUS will likely be larger than 2 mm.
Thank you again for your comments and taking the time to review our work.
Dear reviewer,

thank you for your comments and ideas on how to improve the work.

Evaluating the method on another, public dataset would indeed be beneficial.
The difficulty lies in finding or creating a 4D MRI dataset with sufficient temporal (and spatial) resolution.
Due to the nature of MRI, existing 4D MRI datasets typically use retrospective sorting, potentially in combination with heavy undersampling, to achieve a high enough temporal resolution to visualize different phases of a respiratory cycle.
However, these datasets fail to capture the variability between respiratory cycles and, more crucially, simulating our method on them requires a higher temporal resolution still, so heavy interpolation between the already low-quality images is necessary, lowering the quality and realism of the data further.
Therefore, we have decided not to include such a benchmark dataset.
However, we agree that it would be a valuable bonus to apply our method to a suitable benchmark dataset in the future.
Thank you again for your time and feedback.
Dear reviewer,

thank you very much for your thorough review and detailed feedback!

The pixel intensities for a given organ vary between sequences based on the observed intensity distributions.
However, each pixel within an organ was assigned the same intensity value.
This is a simplification introduced by the XCAT phantom and one of the reasons we add noise during the data generation.
The noise is added on the k-space data but this also results in noise in the image domain.
This noise is global and not sampled from tissue-specific distributions, but it still adds variety within an organ.

The spherical shape and therefore missing rotations are an unfortunate limitation of the XCAT phantom.
Typical rotations of tumors in the liver are small, but important to estimate, especially for large, irregularly shaped tumors.
In our previous work, a separate network was successfully used to estimate the rotation and since rotations here would mainly be induced by the periodic respiration, we are confident that the RNN could learn this relationship in 3D as well.
However, we are aware of this limitation and are considering ways to simulate more realistic data in the future.

The number of sequences per fold in section 2.3 was indeed a mistake and should have been 20.

Ablation studies for the spoke halves and subsequences are a good idea.
Hence, we have performed these experiments and added the results to the appendix.
Using half spokes as input leads to similar results for a TR of 40 but improves the estimations for a TR of 4, compared to using the full spokes.
The subsequence learning has a much larger impact.
Without it, the RNN does not learn much at all for a TR of 4.
Again, the additional variation with a TR of 40 leads to less of an impact of this, but it is still significant.

We agree, that the analysis of the results was not as thorough as it should have been.
We have added additional plots and analysis to the appendix to better understand the observed errors.
This includes the suggested box plots to visualize the sequence-to-sequence variation, as well as a plot of the error over the respiratory phase.
The box plots show that the MAE varies greatly between sequences, but that the RNN and the hypothetical perfect estimator perform similar on the same sequence for a TR of 4.
However, the respiratory phase affects both methods very differently and is another source of variation in the MAE.
While the RNN makes the largest errors around the end full inhale phase when the motion amplitude is the largest, the perfect estimator makes the largest errors between full inhale and full exhale phases, where the speed of motion is fastest.

We have also clarified the comparison to the perfect estimator in the text a bit to highlight what its results mean and what the advantage of the low latency of our method is.

Once again, thank you for your helpful comments. Our manuscript has greatly benefited from them.