Motion correction of free-breathing magnetic resonance renography using model-driven registration

Introduction Model-driven registration (MDR) is a general approach to remove patient motion in quantitative imaging. In this study, we investigate whether MDR can effectively correct the motion in free-breathing MR renography (MRR). Materials and methods MDR was generalised to linear tracer-kinetic models and implemented using 2D or 3D free-form deformations (FFD) with multi-resolution and gradient descent optimization. MDR was evaluated using a kidney-mimicking digital reference object (DRO) and free-breathing patient data acquired at high temporal resolution in multi-slice 2D (5 patients) and 3D acquisitions (8 patients). Registration accuracy was assessed using comparison to ground truth DRO, calculating the Hausdorff distance (HD) between ground truth masks with segmentations and visual evaluation of dynamic images, signal-time courses and parametric maps (all data). Results DRO data showed that the bias and precision of parameter maps after MDR are indistinguishable from motion-free data. MDR led to reduction in HD (HDunregistered = 9.98 ± 9.76, HDregistered = 1.63 ± 0.49). Visual inspection showed that MDR effectively removed motion effects in the dynamic data, leading to a clear improvement in anatomical delineation on parametric maps and a reduction in motion-induced oscillations on signal-time courses. Discussion MDR provides effective motion correction of MRR in synthetic and patient data. Future work is needed to compare the performance against other more established methods.


Introduction
Quantitative imaging biomarkers are measured by acquiring images with different acquisition parameters, and fitting a physical model of the ensuing signal changes to these data. For body areas like the abdomen and thorax, motion correction is often necessary to avoid significant artefacts in the parametric maps. Free-breathing magnetic resonance renography (MRR) [1], or dynamic contrast-enhanced MRI (DCE-MRI) of the kidney, exhibits a combination of characteristics that present extremely challenging conditions for motion correction. This includes major reversals in image contrast and large motion amplitudes compared to the internal structures of the target organ. As a result, and despite intensive research in the area, a sufficiently robust motion-correction approach for MRR has not yet been identified [2][3][4][5][6].
Standard pairwise image co-registrations between the images in the series [7,8] or with a single reference image 1 3 is often ineffective in quantitative imaging due to the large changes in image contrast. The problem can be solved by adopting a more group wise perspective and considering the motion correction of all images in the series as a joint optimization problem. Two different classes of methods can be distinguished: (1) Model-driven registration (MDR) methods generate target images by fitting the physical signal model to the data [9] and perform a (usually pairwise) coregistration between source images and targets. (2) Modelfree registration methods do not rely on the physical signal model and instead use more general signal-processing concepts to provide a cost function for the groupwise co-registration problem [10][11][12][13][14][15][16][17].
The attraction of model-free registration lies in the fact that the same algorithm can be applied to any quantitative imaging method, and the results are not affected by biases in the physical signal model [18]. On the other hand, motion correction is not an aim in itself in quantitative imaging, but part of a processing pipeline that will ultimately always apply a signal model. Hence in this context, model-free registration can impose additional constraints that may bias the solution unnecessarily. For instance, the assumption that breathing does not affect the principal components of a signal is not necessarily true. In particular, this may be the case when data are acquired rapidly in free-breathing and motion creates a large periodic oscillation on the signal.
MDR does not impose additional assumptions beyond those implicit in the physical signal model, and may therefore be a more suitable motion-correction approach in quantitative imaging [19]. Evidence shows that MDR performs well in a wide range of application areas, but the majority of these involve data with relatively small amounts of motion, such as breast [9] and brain [20][21][22] (immobilised by dedicated coils), ECG-triggered cardiac or abdominal data in breath hold [23][24][25][26][27][28][29] or lower abdomen [30][31][32] (limited bulk motion and peristalsis suppressed with antispasmodics [33]). Only two studies so far have applied MDR in free-breathing [19,34], both on abdominal diffusion-weighted MRI which does not exhibit the large reversals in contrast shown by T1-mapping or DCE-MRI. MDR has not yet been evaluated as a motion correction method for MRR.
The aim of this study was to investigate whether MDR is a suitable approach to correct for breathing motion in high temporal resolution MRR acquired under free-breathing. The standard formulation of MDR was generalised to cover this type of problem, and MDR was evaluated on a kidney-mimicking digital reference object (DRO), as well as on patient data using multi-slice 2D and 3D acquisitions. The algorithm, DRO and patient data are freely available as supplementary material to allow independent verification of the conclusions and serve as benchmarks for future developments in motion correction (github.com/plaresmedima/ Flouri et al. 2020).

Theory
This section introduces basic theory and notations, and generalises conventional MDR for application to MRR with linearised kinetic models. The theory is introduced in the most general context possible to allow transfer of experience to other applications of quantitative imaging and DCE-MRI.

DCE-MRI signal model
The basic approach to modelling of MRR is well-known [1] and follows standard principles of DCE-MRI, but some of the usual distinctions are irrelevant for the purposes of motion correction. We repeat here the main definitions to unambiguously fix notations and assumptions.
DCE-MRI measures S(x, t) in location x as a function of time t after injection of a contrast agent. For the purposes of this study, we will assume that tissue concentrations are small enough so that S(x, t) is approximately linear in the concentration [35]: Here, (x) = r 1 S 0 (x)T 10 (x) depends on the relaxivity r 1 of the contrast agent, the pre-contrast relaxation time T 10 and on the pre-contrast signal S 0 . We will also assume that the tracer distributes over at most two distinct compartments with a global input function c a (t) . In that case, C(x, t) is a convolution of c a (t) with a bi-exponential: Here, F i (x) and T i (x) are four unknown parameter fields. The model is easily generalised to more complex systems by adding terms in the sum, one for each compartment [21,36]. Inserting Eq. (2) into Eq. (1) provides the main model equation for two-compartmental DCE-MRI: Here, we absorbed the calibration map (x) into a re-definition F i → F i , as the parameter interpretation is not relevant for motion correction purposes.
To avoid a dependence on initial values and speed up model fitting, a linear formulation of the two-compartment model can also be used [37]: Here, the superscripts '(1)' and '(2)' refer to single and double integration over time, respectively. The unknown model parameters are the four fields i (x) and i (x) , which are directly related to the four parameters F i (x) and T i (x) [37]. Inserting Eq. (4) into Eq. (1) leads to the signal model for linear DCE-MRI: As above we absorbed (x) into a re-definition i → i .

Free-form deformation model (FFD)
Consider a set of images S(x, t) acquired at different times t and corrupted by motion. If motion is modelled by a coordinate transformation x → D(x, t) at each time point t (the deformation field at t ), then motion-free images S D (x, t) can be derived as: We are free to define which is considered the motionfree state D x, t f = x . In a free-breathing protocol, it is convenient to choose t f in end-expiration as this typically covers the largest part of a breathing cycle.
The breathing motion in this study is modelled with free-form deformations (FFD) [7]. The deformation field D(x, t) is defined at any by interpolating vectors D j (t) = D x j , t on a rectangular grid of control points x j with grid spacing Δ = Δ x , Δ y , Δ z : In this study, we use linear interpolation to speed up the computations, in which case the weighting functions W Δ j (x) = W Δ (x − x j ) are defined as: with w(u) = 1 − |u| for |u| ≤ 1 and w(u) = 0 otherwise. The coordinate transformation (Eq. 6) is applied by interpolating the images between the values S l (t) = S(x l , t) at voxel centres x l . If the voxel dimensions are δ, the signal at any location is given by: The deformed image (Eq. 6) can then be parametrized in terms of D j (t) by substituting x → D(x, t) in the righthand side of Eq. (9) and using Eq. (7):

Standard MDR (fixed target)
Consider an arbitrary quantitative imaging method where signals S(x, t) are measured at locations and at times t with different signal parameters p(t) . Depending on modality, p can represent time points directly (e.g., DCE-MRI or Dynamic PET), b-values and gradient directions (diffusion MRI), inversion or echo times (MRI relaxometry), phaseencoding directions (4D Flow), etc. Standard MDR assumes there exists a signal model ∑ (P(x), p) that describes motion-free images in terms of N unknown parameter maps P(x) = P 1 (x), … , P N (x) : The left-hand side is here the motion-corrected image as defined in Eq. (6). Equation (11) is the defining equation of standard MDR, providing a direct link between deformation fields D(x, t) and parameter maps P(x). Equation 11 is solved by initialising D(x, t) = x and iterating the following two steps until convergence: (1) solve for P at fixed D; (2) solve for D at fixed P. This breaks up the problem in a series of smaller problems that can be addressed with existing methods: solving for P at fixed D amounts to fitting the signal model to a series of images S D (x, t) ; solving for D at fixed P presents a standard pairwise image co-registration problem for each t that can be solved with common software packages.
While image co-registration generally requires normalised cost functions to account for differences in contrast and intensity between source and target, in MDR, a least squares cost function can be used because the source images S D (x, t) and target images Σ(P(x), p(t)) have similar intensity and contrast. Since the physical signal models are typically also solved in the least squares sense, this allows us to formulate MDR as a joint optimization problem: Since both steps of the iteration then optimize the same cost function, convergence is guaranteed. The use of a least squares cost function also offers a significant computational advantage because it affords an analytical expression for the partial derivatives G j (t) of the cost function with respect to each D j (t) [38]: Here (∇S) D is the FFD of the image gradient ∇S , R(x, t) is the residual S D (x, t) − Σ(P(x), p(t)) , and N j is the support of the wavelet W Δ j (x) (Eq. 8). Equation (13) can be computed efficiently. The weights W j (x) only depend on the resolution level Δ and can be pre-computed for each control point j . The product R(x, t) (∇S) D (x, t) needs to be computed only once for all control points j , and the summation over x can be restricted to a small neighbourhood N j around each control point j (Eq. 8).

Generalised MDR (moving target)
A hidden assumption in the formulation of standard MDR is that the signal model Σ(P(x), p) is not a function of the measured signals S(x, p) themselves. In general, however, this is not necessarily the case.
Consider the example of the DCE-MRI model in Eq. (3). If S 0 (x) is treated as an unknown model parameter, then Eq.
, and the assumptions of standard MDR are fulfilled. In practice, however, nearly all current implementations of DCE-MRI derive S 0 (x) from the data by averaging the pre-contrast signals S(x, 1), … , S x, n 0 . In that case, , t becomes a function of the signal at different times. This is also true when a linear implementation is used as in Eq. (5), which depends explicitly on S(x, t) = even if S 0 (x) is treated as a free parameter.
The implications are significant. If the signal model depends on the signal itself, then Eq. (11) must be recast as follows 1 : Solving this equation for D at fixed P no longer presents a standard image co-registration problem because the righthand side also depends on D-presenting effectively a moving target. An additional problem is that a deformation at any single time t can now affect the equation at any other time t ′ . This implies that co-registrations at different times can no longer be treated as independent problems, and the problem requires a groupwise optimization of the cost function: As shown in the appendix, the gradient of the cost function (Eq. 15) has the same form as the special case (Eq. 13), after substituting a generalised residual R(x, t) → R(x, t)R: 2 We used the notation t R (x, s) for the partial derivative of The need for a groupwise optimization in generalized MDR comes with a significant computational overhead. However, in the case of DCE-MRI with a linear model (Eq. 5), the assumption t R (x, s) ≈ ts may well present a reasonable approximation to the exact gradient. Defining a matrix R(x) with elements ( R) ts (x) = t R (x, s) , and writing I for the identity matrix, ΔtMS for numerical integration of S over time, and M 0 S for the K-element vector S 0 , … , S 0 T , we can write the partial derivative of the residual as: where we used the definition of i and i [37] to define: In DCE-MRI, the time step ∆t is always chosen to be significantly shorter than the mean transit times of the system, so that Δt ≪ T 1 , T 2 and therefore ∈ 1 , ∈ 2 are small quantities. Since further, the non-zero elements of the matrix M 0 are 1/n 0 (with n 0 ∼ 10-20 the number of pre-contrast images), the approximation R ≈ I may be justified. In that case, the dependency on D of the right-hand side in Eq. (14) can be ignored and the equation can be solved in the same way as standard MDR. This approach has been adopted in the current study.

Implementation
MDR was implemented with the linear DCE-MRI signal model (Eq. 4). A multi-resolution strategy was adopted where the deformation fields D j (t) are initially determined over a coarse grid of control points with spacing ∆ equal to the field of view. The solutions D j (t) are then interpolated to initialise deformation vectors at a finer grid ∆ ⟶ ∆/2. This process is iterated until a user-defined minimum grid spacing ∆ min . Deformation fields at the lowest resolution were initialised to no deformation: At each fixed resolution level ∆, the interpolation weights (Eq. 8) were pre-computed and the following two steps were iterated: (1) pixel-by-pixel linear least squares fitting of the four model parameters i (x), i (x) to the deformed images S D (x, t) [37]; (2) time-by-time fitting of the deformation parameters D j (t) using gradient descent with a back-tracking line search and gradients calculated with Eq. (13). The line searches were preconditioned by using the step size of the previous line search to initialise the next. The gradient descent was stopped if the line search returned a maximum change in D j (t) less than a user-defined tolerance δ min . The iterations at a given resolution level were stopped if none of the time points produced a change larger than δ min . Optimal values for the minimum grid spacing ∆ min and the tolerance δ min were found for a test case for each type of data in two steps. First, the grid spacing ∆ min was varied at the lowest tolerance considered (δ min = 0.1 pixels) and an optimal grid spacing was identified. Second, to minimize computation times, the tolerance δ min was varied between 1.0 and 0.1 for the optimal grid spacing ∆ min selected in the previous step. The smallest tolerance beyond which no substantial improvement was apparent was identified.
MDR was implemented in IDL 6.4 (Exelis VIS, Boulder, CO) and run on a standard laptop PC with a 2.7 GHz Intel Core processor and 16 GB memory.

Data and simulations
MDR was evaluated using a combination of 2D simulated data and patient data acquired in 2D and 3D. All patients provided written informed consent and the studies were approved by the local research ethics committees. The 3D data were co-registered with 3D MDR where D j (t) has 3 independent components. 2D data were co-registered with 2D MDR.

Digital reference object (DRO)
A 2D DRO with kidney-like structures was used to generate motion-corrupted MRR images (120 dynamics, 1.1 s intervals, matrix size 135 × 135). The "kidneys" were modelled as concentric ellipses with different values for the parameters F P and T T to model cortical, medullary, and pelvic structures. The parameter values in the different regions range from 30 to 300 mL/min/100 mL (plasma flow F P ), 6-20 s (plasma mean transit time T p ), 40 to 60 mL/min/100 mL (tubular flow F T ) and 90 to 300 s (tubular mean transit time T T ). A literature based c a (t) was used [39], pre-padded with zeroes to create a 15 s baseline. Data were generated with Eq. (3) and S 0 (x) = 1 at pseudo-continuous temporal resolution 0.1 s, and then down-sampled to 1.1 s.
Three different types of motion were applied: no motion, rigid motion through sinusoidal vertical shifts with amplitude 12 pixels and period 4 s, and non-rigid motion derived from the deformation fields measured on a 2D clinical data set. The motion amplitudes were scaled up compared to the clinical data to test the performance of the approach under challenging conditions. The no-motion case is included to present a best-case scenario for MDR. Gaussian noise was added to the signal with varying standard deviations (SD) to test the noise sensitivity. Noise levels are expressed in terms of the contrast-to-noise ratio (CNR) defined as max (c a )/SD.
A particular issue in the design of DRO's for motion correction is the choice of the ground truth. A result where the motion is frozen in a position that is different from that of the "ground truth" parameter maps is not necessarily incorrect as long as the frozen position is a natural breathing state. The problem is addressed in the DRO by applying the deformation fields to the "ground truth" parameter maps as well, and using as actual ground truth for error quantification the breathing state that minimises the difference between the reconstructed plasma flow and moving plasma flow map.

2D clinical data
2D MRR data were taken from five consecutive subjects enrolled in an ongoing study into MRI biomarkers of renal fibrosis. The study was approved by the local research ethics committee, and all subjects gave written informed consent (Newcastle and North Tyneside 1 Research Ethics Committee, REC reference 14/NE/1120). MRR was performed in free-breathing using a 3.0 T scanner (Philips Achieva, Best, The Netherlands) and a 2D saturation-recovery turbo-flash sequence with linear encoding of k-space. A 2-channel body transmit coil was employed for homogeneous signal transmission and data were acquired using an 18-channel torso receive coil. Four slices (3 coronal, 1 axial) were acquired at a temporal resolution of 1.1 s. Other imaging parameters were as follows: acquisition matrix 116 (phase) × 135 (read), number of dynamics 250, echo time 1.63 ms, repetition time 3.6 ms, bandwidth 900 Hz, saturation recovery time 148 ms, flip angle 12 • , slice thickness 7 mm, sensitivity encoding (SENSE) factor 2.4, field of view 375 × 440 mm 2 , in-plane resolution 3.2 × 3.2 mm 2 , reconstructed matrix 480 × 480. The contrast agent Gd-DOTA (Dotarem, Guerbet, France) was injected with a half dosage of 0.05 mmol/kg body weight. The input c a (t) was measured in the aorta on the axial slice.

3D clinical data
3D MRR data were taken from eight consecutive cases from a study on atherosclerotic renovascular disease [40]. The study was approved by the local research ethics committee, and all subjects gave written informed consent (Oldham Local Research Ethics Committee, REC reference 07/Q1405/21). MRR was performed using a 3.0 T whole body scanner (Philips Achieva, Philips Medical Systems) with a phased-array body coil for signal reception. Subjects were imaged using 3D spoiled gradient echo sequence in the oblique coronal plane. The following parameters were used for the acquisition: repetition time 5.0 ms, echo time 0.9 ms, field of view 400 × 400 × 100 mm, voxel size 3.13 × 3.13 × 4.0 mm 3 , flip angle 17 • , SENSE factor 2, acquisition matrix 128 × 84 × 10, reconstructed matrix 128 × 128 × 20. This resulted in a temporal resolution of 2.1 s/volume. Subjects were given a quarter dose of 0.025 mmol/kg GdDOTA (Dotarem, Guerbet, France) at a rate of 3 ml/s. In all cases, the acquisition and contrast agent injection were initiated simultaneously. The input c a (t) was measured in the aorta between the bifurcations of renal and iliac arteries.

Assessment of registration
Registration quality in the DRO was measured by the percent error in any given field P( ) , quantifying the difference between the reconstructed P r (x) and the ground truth P(x): Bias of the result was quantified as the median E P (x) across all x, and the precision as the width of the 90% confidence interval. In addition to the quantitative metrics, registration quality was also assessed visually by comparing reconstructed parameter maps against exact maps P(x).
The registration quality in the patient data was assessed qualitatively by visual comparison of corrected and uncorrected parameter maps, time curves and model fits in motion-sensitive areas, and time-cut images.
Quantitative assessment of MDR was performed based on manual regions of interest (ROI) corresponding to the left kidney and right kidney. To create the ground truth, the left and right kidney were manually segmented on a single slice of 2D renal data. This was then propagated across all the time frames. The segmentation of the unregistered data was obtained by manually adjusting the position of the ROIs in every time frame to best follow the feature of interest. The same procedure applied to that of unregistered data was followed to obtain segmentation of registered data. Each segmentation was evaluated against the ground truth by calculating the Hausdorff distance using the EvaluateSegmentation tool [41]. A 2D slice was extracted from the middle of the 3D renal data and the same procedure applied to that of 2D renal data was followed to obtain the segmentations of 3D data.
To assess the impact of model bias, MDR was also applied with a simplified 3-parameter modified Tofts kinetic model [42]. The model is known to provide a poor fit to the data in the first pass, which introduces bias in the motion correction.

Statistical analysis
The statistical analysis was performed using RStudio (version 4.0.3, 2020). Normality was assessed with Shapiro-Wilk test. Paired t test was performed to compare the Hausdorff distance between the unregistered and registered data. The 2D and 3D data were grouped together, and the analysis was carried out on the combined data. The significance level was set at 5%. Figure 1 illustrates the source data used in the study. The 3D data showed significant intra-frame artefacts, presenting a particular challenge for motion-correction. Figure 2 shows the effect of varying the deformation field gridsize ∆ min on the parameter maps. For the patient data, the F P maps show a gradual sharpening of the image until the smallest gridsize, but F T maps show that at the smallest gridsizes of 4-and 8-pixel small-scale deformations are induced that degrade the expected anatomical structure. Inspection of the dynamics suggests that at this stage, the deformation field attempts to correct for model errors. The optimal grid size for these data was chosen to 32 pixels (dashed lines) as no obvious further improvements are seen at 16. For the DRO, results improve when the grid size is reduced to 32 but no further improvements are seen beyond that. The effect of model bias does not play a role in the DRO because forward and inverse models were the same. Figure 3 shows the effect of varying the tolerance δ min . The patient data show that the results gradually converge to the result at the lowest tolerance of 0.1, at the cost of significant extra calculation time. With tolerances lower than 0.3, the calculation time increases dramatically from 67 min at 0.3 pixels to 88 min at 0.2 pixels and 489 min at 0.1 pixels. As there is no obvious improvement visible at tolerances lower than 0.3, this was chosen to be the optimal tolerance for these data-saving around 6 h of computation time per slice compared to the lowest tolerance of 0.1. The DRO shows a similar trend, but due to the sharp tissue interfaces, improvements continue to be apparent until a tolerance of 0.2 pixels. Figure 4 illustrates the results of MDR on the simulated data and a single clinical example. Comparison to the uncorrected data in Fig. 1 shows that the algorithm removes the motion effectively, including the exaggerated motion in the DRO. Figure 5 illustrates the effect of motion correction on temporal profiles, model fits and parameter maps for simulated and patient data. In all cases, the time curves show strong breathing-induced oscillations that are substantially reduced after motion correction, leading to a much-improved fit of the model to the data. In simulated and 2D data, the model describes the motion-corrected data almost exactly. In the 3D data, some oscillations remain after co-registration, consistent with the effects of within-frame motion artefacts. On the parameter maps, registration significantly reduces the motion-induced blurring that is visible on the uncorrected maps. This leads to clearer organ boundaries and delineation of internal anatomical structures, such as renal cortex and medulla. Figure 6 shows the effect of registration through time cuts in all five 2D datasets. The time cuts of the unregistered data clearly show that significant motion is present which varies in structure and amplitude between subjects. The registered data demonstrate that the motion has effectively been removed, without affecting the changes in contrast. Typical calculation times for a single 2D slice were around 40 min on a laptop PC. Figures 7 and 8 show the effect of registration through time cuts in all eight 3D datasets. Motion correction in the 3D data is more challenging due to the appearance of significant within-frame motion artefacts, which cannot be removed by image co-registration. Nevertheless, the time cuts clearly indicate effective reduction of motion effects. Visual inspection of the dynamic data demonstrates that kidney motion is effectively frozen, and any residual breathingrelated oscillations are due to the within-frame artefacts. For all 3D data, the grid size was set to 4 pixels and tolerance to 0.2, leading to a calculation time of about 14 h on a laptop PC. Figure 9 shows the effect of registration through binary segmentation masks. Visual inspection of the ground truth segmentation mask demonstrates that kidney motion is efficiently removed. Table 1 presents the Hausdorff distance values used to assess segmentation quality before and after motion correction. Analysis of the 13 patients showed that the Hausdorff distance (HD) was significantly lower in segmentations after motion correction compared to segmentation masks before motion correction (HD unregistered = 9.98 ± 9.76, HD registered = 1.63 ± 0.49, P value = 0.00027). Figure 10 shows the error distribution as measured on the DRO before and after MDR, with CNRs ranging from high to low noise levels. The results in the motion-free case show that MDR improves the precision (width of 90% CI) at high and medium noise levels due to the spatial smoothing effect of the coordinate transformation. Comparing the results without MDR in the presence of rigid and non-rigid motion to the motion-free case, the figure shows that the motion induces additional error-mainly in the form of reduced precision under low noise conditions. At high noise conditions, the error is fully dominated by noise with no obvious added contribution of motion. After MDR, the error in data with rigid and non-rigid motion is reduced and becomes indistinguishable from that measured in motion-free data. Figure 11 illustrates the effect of bias in the signal model. The modified Tofts model assumes infinite flow and therefore provides a poor fit to the earlier time points during the first pass of the contrast agent. As a result, the target for co-registration is a poor match to the actual image in the first pass, resulting in co-registration errors. This is illustrated in the figure: at the times corresponding to the first pass ( t 20 , t 24 , t 25 in the simulated data and t 37 − t 39 in patient data), the co-registered images are deformed.

Discussion
The aim of this study was to investigate if MDR was suitable for motion correction in the challenging application area of free-breathing MRR at high temporal resolution.
A key difference is that previous MDR methods [19,23,25,28,29] use pre-packaged co-registration modules, which does not allow to exploit the inherent symmetries of the problem to eliminate duplicate calculation steps. The MDR implementation proposed in this study uses a tailored FFD-based registration model which enables Fig. 2 Selection of an optimal grid size ∆ min for the example of 2D patient data (top rows) and DRO with non-rigid motion (bottom rows). All calculations are performed at the lowest tolerance of δ min = 0.1 pixels. For each dataset the figure shows the results for plasma flow (F P ) and tubular flow (F T ) at different grid sizes from lowest (left, 4 pixels and 2 pixels, respectively) to highest (right, 256 pixels and 128 pixels, respectively). The unregistered result is shown in the right most column. For the DRO the ground truth is given on the left. Motion corrections are performed over the entire image, but results are cropped to the right kidney to show the detail. The dashed lines show the selected grid size for each case. Dynamics, results for the other two parameters and error maps for the DRO are shown in the supplementary material significant acceleration, allowing a co-registration in 3D with acceptable computation times on a standard PC without dedicated CPU implementations. Results in synthetic and patient data have shown that MDR effectively removes the motion in calculation times that are feasible on standard laptop computers.
MDR suffers from bias in cases where the chosen signal model does not adequately describe the data. Our results have confirmed this effect, showing significant unphysical deformations when an oversimplified model is used. The motion correction, when performed at small enough grid spacing, will in those conditions attempt to correct for errors in the signal model by collapsing or deforming areas of the image where the signal model does not fit the data. In qualitative imaging applications, MDR may therefore not be a preferred approach to motion correction [18]. However, in quantitative imaging a model-error will always bias the results, and this problem is not resolved by eliminating it from the motion-correction process. In fact, as our results show, the ensuing unphysical deformations help to flag up the presence of model error, where this may not be immediately obvious from a goodness-of-fit.
In a recent study, Coll-Font et al. [43] argued that MDR is unsuitable in the context of MRR because the signal models are valid only for the kidney, and that this therefore necessitates a region-of-interest-based co-registration. The predicted issues were not encountered in this study despite the approach of co-registering over the entire field of view. Indeed, what matters in MDR is the fit of the model to the data, and from this perspective, the renal compartmental model is no different than DCE-MRI models for other tissue types in the abdomen. The impulse response function is modelled as a bi-exponential, and it is only the interpretation of the model parameters that is specific to the kidney (Eq. 2). It should be noted though that this may not be true in other application areas. For instance, in cardiac DCE-MRI where an input function is taken in the left ventricle, the model will not be able to fit the lungs, pulmonary arteries or right ventricle as the bolus arrives in those tissues before the left ventricle. In this case, the validity of the model may indeed Fig. 3 Selection of an optimal tolerance δ min for the example of 2D patient data (top rows) and DRO with non-rigid motion (bottom rows). All calculations are performed at the optimal grid size of 32 pixels (see Fig. 2). For each dataset, the figure shows the results for plasma flow (F P ) and tubular flow (F T ) at different tolerances from lowest (left, 0.1 pixels) to highest (right, 1 pixel). For the patient data, the calculation times for each tolerance level are also shown. Motion corrections are performed over the entire image, but results are cropped to the right kidney to show the detail. The dashed lines show the selected tolerance for each case. Dynamics, results for the other two parameters and error maps for the DRO are shown in the supplementary material be more limiting, and the application may benefit from a model-free registration followed by a region-of-interestbased model fitting. An alternative solution though may be the use of MDR but with an input function selected in the pulmonary artery.
We used a FFD model for the breathing motion, as this is very generally applicable, can be easily adapted to motion at different scales by changing the smallest grid spacing, and is comparatively efficient computationally. However, other motion models may be used and could offer a benefit. For instance, a number of recent studies using MDR applied diffeomorphic transformations [30,32], which are also free-form but derive deformations from stationary velocity fields. This ensures that the ensuing deformations are invertible, which is a necessary condition of any physical deformations. Diffeomorphic registrations will help to avoid unphysical deformations, producing more physical results even in the presence of model bias.
This may be suitable in scenarios where an accurate physical model is not available or would be underdetermined by the available data.
In this study, we chose to use a newly built implementation of MDR without using pre-packaged solutions. For prototyping, this has the advantage of transparency, allowing a dissection of the different components, identification and precomputation of repeated calculation steps, and easy evaluation of alternative algorithm architectures. However, this also implies there is significant scope for acceleration using well-known recipes, such as stochastic gradients, preconditioner estimations for the line search [44], or optimizing at reduced resolution. Co-registrations at different times can be parallellized, and well-optimized CPU-type implementations of FFD can be used [45]. More generic registration packages, such as Elastix [46] and ANTs, [47] can easily be integrated to evaluate alternative approaches. Machine learning may also lead to step changes in computation times [48]  The plots show the signal-time curves (dashed line) and the model fit (solid line) before and after motion correction. The plasma flow (F P ) maps before and after motion correction are also presented and can presumably be implemented using unsupervised neural networks [49].
A separate potential route for improvement could involve a groupwise approach with the exact Eq. (16) rather than the approximation ( R) ts (x, s) = ts . This was explored in preliminary experiments and abandoned after these demonstrated extremely slow convergence without obvious benefit. Nevertheless, the approach has not been evaluated systematically and can potentially be improved using strategies, such as joint optimization of deformation fields and parameter maps. The alternative approach of using the analytical model solution (Eq. 3), which fulfills the condition ( R) ts ( , s) = ts , was also explored. The approach was abandoned after preliminary experiments demonstrated salt-and-pepper noise that dominated the cost function in the co-registration step. However, potential solutions, such as improved non-linear optimization or de-noising, have not been explored.
While the focus of this study is on retrospective motion correction by image registration, this is not the only means   of correcting or reducing motion in medical imaging. Other approaches include motion-corrected compressed sensing [50,51], prospective or retrospective gating, breath holding, using k-space trajectories that are less motion-sensitive [52], or a combination of the above [17]. These methods can be combined with MDR, and there may also be a benefit in integrating MDR with image reconstruction to produce a joint model-driven reconstruction [53,54] and motion correction. Finally, while we have demonstrated that MDR effectively corrects for breathing motion in MRR, future studies are needed to determine whether MDR improves on alternative model-free approaches that have been proposed in this context [15,43].

Conclusion
MDR provides effective motion correction in the challenging application of free-breathing MRR at high temporal resolution. Future studies are needed to determine whether MDR improves the results compared to alternative co-registration methods proposed in this context.  10 Error distribution before and after motion correction for the simulated data at CNR from 10 0 (high noise level) to 10 4 (low noise level). The columns show the two parameters and the rows show different motion types: no motion (top row), rigid motion (middle row), non-rigid motion (bottom row). The circles indicate the median relative parameter error for all the pixels in the image, and the error bars represent the 90% confidence intervals. For reasons of clarity, only two of the parameters are displayed; the trends in the other two parameters were similar