Correcting inter-scan motion artifacts in quantitative R1 mapping at 7T.

PURPOSE: Inter-scan motion is a substantial source of error in R1 estimation methods based on multiple volumes, for example, variable flip angle (VFA), and can be expected to increase at 7T where B1 fields are more inhomogeneous. The established correction scheme does not translate to 7T since it requires a body coil reference. Here we introduce two alternatives that outperform the established method. Since they compute relative sensitivities they do not require body coil images. THEORY: The proposed methods use coil-combined magnitude images to obtain the relative coil sensitivities. The first method efficiently computes the relative sensitivities via a simple ratio; the second by fitting a more sophisticated generative model.
METHODS: R1 maps were computed using the VFA approach. Multiple datasets were acquired at 3T and 7T, with and without motion between the acquisition of the VFA volumes. R1 maps were constructed without correction, with the proposed corrections, and (at 3T) with the previously established correction scheme. The effect of the greater inhomogeneity in the transmit field at 7T was also explored by acquiring B1+ maps at each position.
RESULTS: At 3T, the proposed methods outperform the baseline method. Inter-scan motion artifacts were also reduced at 7T. However, at 7T reproducibility only converged on that of the no motion condition if position-specific transmit field effects were also incorporated.
CONCLUSION: The proposed methods simplify inter-scan motion correction of R1 maps and are applicable at both 3T and 7T, where a body coil is typically not available. The open-source code for all methods is made publicly available.

INTRODUCTION

Quantitative MRI, and the push toward in vivo histology, aim to extract tissue‐specific parameters from a series of weighted volumes. For example, the longitudinal relaxation rate, , which is sensitive to important biological features, such as myelin and iron content, can be quantified with the variable flip angle (VFA) approach, for example, References 2, 3. A common assumption when computing quantitative metrics is that certain multiplicative factors, such as the signal intensity modulation imposed by the receiver coil's net sensitivity profile, are constant across the weighted volumes. However, this is invalid if motion occurs between the volume acquisitions. In the case of neuroimaging, rigid body co‐registration can be used to realign the brain but will not correct for the differential coil sensitivity modulation, which in maps computed with the VFA approach can lead to mean absolute error approaching 20%.

A correction scheme has previously been proposed by Papp et al. and validated for mapping at 3T. The position‐specific net receive sensitivity is estimated from two rapid low‐resolution magnitude images, received on the body and array coils respectively prior to each VFA acquisition. The more homogeneous profile of the body coil is used as a reference to compute the net receiver sensitivity, which is then removed from the VFA acquisitions. This approach effectively assumes that the body coil's modulation is consistent across volumes instead of that of the array coil. This in itself is a potential limitation, as is the general unavailability at body coils at higher field strengths.

Here we propose an alternative whereby we estimate the relative sensitivity between volumes. This approach does not fully remove the receiver's sensitivity modulation but does remove the bias that differential modulation introduces in quantitative metrics. Only the calibration images obtained with the array coil are required, that is, less data than the originally proposed method. To validate the approach, we focus on maps computed with the multiparameter mapping (MPM) protocol. We first compare performance with the established method of Papp et al. at 3T and then demonstrates a reduction of inter‐scan motion artifacts at 7T under a range of different motion conditions. We further demonstrate that, unlike at 3T, the transmit field also exhibits substantial position‐specific variability at 7T. As a result, the most precise estimates were obtained by accounting for both position‐specific transmit and receive sensitivity effects.

While we validate this approach in the context of mapping, it has much more general potential and can be applied to other mapping methods that combine data from multiple volumes.

METHODS

Theory

mapping can be achieved by acquiring spoiled gradient echo volumes with at least two different flip angles. , At a given spatial location, the image intensity, I, for a given nominal flip angle is: where s is the receive sensitivity, is the proton density, is the transmit field, is the longitudinal relaxation rate, TR is the repetition time, and k indexes the VFA acquisition. Co‐registration allows for inter‐scan motion by realigning anatomical structure across acquisitions. Under the small flip angle approximation, with two nominal flip angles (k = {1,2}), can be computed as follows: Typically, it is assumed that and the sensitivities simplify out. However, this assumption is invalid if inter‐scan motion has occurred leading to substantial bias in estimates. This can be avoided by accounting for the relative sensitivity across positions: . Substitution for in Equation (2) gives: The method of Papp et al. did not include the relative sensitivity but referenced to an additional calibration image acquired on the body coil, assuming that the body coil modulation was position‐independent.

It is commonly assumed that the transmit field is sufficiently smooth as to be considered position‐independent, that is, , such that: However, in this work we show that this assumption does not hold at 7T, and that incorporating position‐specific transmit field estimates maximises the precision of .

Ratio approach

The calibration data used to correct inter‐scan motion artifacts comprised rapid low‐resolution, coil‐combined magnitude images acquired immediately prior to each high resolution VFA acquisition. These images, , assumed to have been rigidly co‐registered to the same space, can be written as the product of a common image r and a net sensitivity field . The relative sensitivity, can be computed with respect to one of the calibration acquisitions, used as a reference:

Dividing each VFA acquisition by its relative sensitivity results in a common modulation, , which, although less homogeneous than the body coil used by Papp et al., more faithfully restores the validity of assuming common modulation when computing .

The assumption that r is common, such that = holds only if there are no position‐specific transmit field effects. Simulations were used to explore the validity of this assumption.

Generative approach

This ratio approach risks noise amplification, particularly in regions of low signal‐to‐noise ratio (SNR). This is combatted by isotropically smoothing and before taking their ratio. A potentially more robust alternative is to cast the computation of the relative coil sensitivities, and a common image modulated by them, as an inference problem in a probabilistic generative model of that incorporates noise and can also embed knowledge about the spatial smoothness of the sensitivities. This generative modeling approach allows coils with arbitrary sensitivity to be incorporated, for example, coils with more (array) or less (body) spatial variation, or both concurrently (“array + body”) if available. A priori knowledge about the expected smoothness of the sensitivity can be incorporated at the level of coil type (body vs. array) via appropriate tailoring of a regularisation parameter, . Images acquired with the body coil will have a flatter sensitivity field modulation, which can be incorporated by setting . Full details are given in Appendix A.

Experiment

Participants

One participant (female, 31 years) was scanned at 3T (MAGNETOM Prisma, Siemens) using a body coil for transmission and either the body coil or a 32‐channel head array coil for reception. Three additional participants (2 female and 1 male; 32–41 years) were scanned at 7T (MAGNETOM Terra, Siemens) using an eight‐channel transmit, 32‐channel receive head array coil (Nova Medical) in a quadrature‐like (“TrueForm”) mode. All data were acquired with approval from the UCL research ethics committee.

MPM datasets

MPM data were acquired using a multi‐echo spoiled gradient echo sequence with flip angles of 6° (PD‐weighted, “PDw”) and 26° (T1‐weighted, “T1w”), a TR of 19.5 ms and an RF spoiling increment of 117° with a total dephasing gradient moment per TR of 6. Eight echoes were acquired with TE ranging from 2.56 to 15.02 ms in steps of 1.78 ms using a bandwidth of 651 Hz/pixel. Data were acquired with a nominal 1 mm isotropic resolution over a field of view of 160 mm right‐left and 192 mm in the anterior‐posterior and superior‐inferior directions. Elliptical sampling and partial Fourier, with factor 6/8 in each phase‐encoded direction, were used to accelerate the acquisition, leading to a scan time of 5 min per volume. A map was estimated by acquiring a series of spin and stimulated echoes using previously described 3T and 7T protocols. , These data were acquired with 4‐mm isotropic resolution resulting in a total acquisition time of 3 min 48 s, and a further 1 min for mapping.

For inter‐scan motion correction, additional single echo acquisitions were acquired prior to each VFA acquisition to facilitate estimation of the relative receive field across positions. These data were acquired with a flip angle of 6°, TE = 2.4 ms, TR = 6.5 ms, a bandwidth of 488 Hz/pixel and no acceleration schemes. At 3T, these data were acquired, receiving sequentially on the array and body coils, with 8‐mm isotropic resolution leading to a scan time of 6 s per volume. To capture the greater spatial variation in the net sensitivity field at 7T, the resolution was increased to 4‐mm isotropic leading to a scan times of 18 s per volume, but acquired only on the array coil due to the absence of a body coil.

Motion conditions

Two MPM datasets were acquired to define baseline reproducibility. Participants were then instructed to move to a new, arbitrary position within the confines of the coil. A localizer was acquired and the field of view repositioned as necessary to ensure appropriate brain coverage in the new position. A third MPM dataset was then acquired.

per contrast

For each map computed from data across two positions, that is, with inter‐scan motion, two different corrections for transmit field inhomogeneity were performed. The first assumed the transmit field was identical across head positions (Equation 4) and in the same position as the PD‐weighted volume. The second used position‐specific maps (Equation 3).

analysis

maps were computed using the hMRI toolbox, which extrapolates the VFA signals to a TE of 0 ms, uses the small flip angle approximation and corrects for imperfect spoiling. The estimation always used the map acquired in the space of the PDw acquisition. Maps were computed with and without inter‐scan motion using all possible PDw and T1w combinations. To ease comparisons, all maps were constructed in the space of the first PDw volume. Rigid transformations between all volumes (calibration data and VFA acquisitions) and the first PDw volume were estimated using SPM12 (Wellcome Centre for Human Neuroimaging) having first corrected for intensity nonuniformity and skull‐stripped the images. maps were computed with and without the proposed inter‐scan motion correction schemes. At 3T, maps were also computed with the method of Papp et al. An isotropic kernel of 12‐mm full‐width‐at‐half‐maximum—the default in the hMRI toolbox—was used to smooth the calibration images prior to computing the relative (proposed) and absolute (Papp et al.) sensitivities. The generative modeling approaches used , , and 15 iterations.

Error metric

Three (two from position one and the third from position two) maps, with no additional inter‐scan motion corrections applied, were averaged to produce a “ground‐truth” map, . For each participant, all of the available maps across conditions (motion/no motion) were assessed against this reference to quantify the error, and its variability. The set of maps used to compute the reference were also segmented to create a mask selecting those voxels with a mean probability of being in WM, GM, or CSF greater than 50%. For participant 3, the cerebellum was excluded, using the SUIT toolbox in SPM, as a result of mapping failure caused by excessively large off‐resonance. For the N voxels within the resulting participant‐specific mask, the mean absolute error, MAE, for each map was computed with respect to the ‘ground‐truth’ map, as: These errors are reported as percentages.

Simulation study

Validity of assumptions

Equation (5) assumes that the calibration data are insensitive to changes in the transmit field across positions such that . Here we test the validity of this assumption via simulation. Under the small flip angle approximation and allowing for position‐specific transmit and receive fields, the calibration images can be written as: Here k indexes the repetition of the acquisition, that is, the calibration data for each high resolution VFA acquisition, and c denotes the calibration‐specific sequence settings. Considering the ratio method for simplicity, can then be written more fully as: For the ratio of the calibration images, , to equal the relative sensitivity, , we require: We note that the Ernst angle is , and rearrange to give: This condition is met when , that is, there is no change in transmit field, or when . These conditions are highlighted in Figure 1 which shows as a function of for . Deviation of from 1 is within 3% for a broad range of values ( and centred on both acquisitions being at the Ernst angle.

FIGURE 1

Misestimation of the true relative sensitivity () by the ratio of calibration images (), as a function of transmit fields. Colours encode , which is 1 when either (red dotted line) or (white dotted line). Isocontours (black solid lines) of are overlaid

Theoretical error

Numerically, errors in the estimates were computed for , and the empirically observed range of relative transmit and receive fields. The median proportion of error arising from transmit or receive field changes was computed over this 4D parameter space.

RESULTS

Exemplar images, relative sensitivities, and results from the generative modeling are shown in Figure 2. The and error maps obtained at 3T and 7T are shown in Figures 3 and 4, respectively. The means and standard deviations of the MAE are reported in Table 1. The differential impact of correcting for transmit and receive field effects is illustrated in Figure 5. This shows and error maps without motion, and with motion having implemented (i) no correction, (ii) correction only for receive field effects, (iii) only for transmit field effects, or (iv) for both effects in combination.

FIGURE 2

3T Example. The acquired calibration images, have different orientation due to participant movement between acquisitions (1st row). Co‐registration can align the images spatially, but does not correct for their differential sensitivity field modulation (second row), visible via their ratio, . The method from Papp et al. estimates and corrects this modulation using an additional body‐coil image (third row). When one is not available, relative signal differences can be corrected for using the relative modulation (fourth row). Alternatively, the joint log‐likelihood of a generative forward model that embeds the spatial transformation from a mean image, r, to native space can be maximized to determine the mean image and modulating sensitivities, , that best explain the acquired images (fifth row). The generative modeling approach produces a similar relative modulation () but allows for the corrected images to have the minimal modulation of the mean image

FIGURE 3

Results at 3T. The first row shows example maps constructed with each method. The second row shows (normalized) error maps with respect to the ground‐truth map . The third row shows histograms of the ground‐truth (solid line) and corrected (filled area) within the GM (green) and WM (purple); these histograms display probability distributions and therefore integrate to 1. “Generative / only array” used only the array coil images in the generative modeling framework, whereas “Generative / array+body” incorporated both the array and body coil images using coil‐specific regularization for the smoothness of the sensitivity modulation

FIGURE 4

Results at 7T. The first column shows the ground‐truth map for the three participants sorted based on the magnitude of inter‐scan motion. Uncorrected and corrected maps are shown in the middle with the corresponding (normaliszd) error maps with respect to the ground‐truth on the right. Correction is only applied for net receive sensitivity modulation and not for transmit field effects. Rows 1 to 3 of the figure correspond to datasets 1, 2, and 3 as reported in Table 1. The net motion is summarised as the root‐sum‐of‐squares, across the three orthogonal axes, of the translations or rotations independently. In the absence of overt motion, the average displacements between the VFA scans, across the group, were 1 mm and 0.6° for translations and rotations respectively, with a maximum translation of 1.6 mm and a maximum rotation of 1.0°

TABLE 1

MAE (mean s.d. across repeats, in %) with respect to the average reference

	Dataset	Motion	No correction	Ratio	Generative	Generative (array + body)	Papp et al.
3T	#1	No	3.0 ± 0.1	3.0 ± 0.1	3.0 ± 0.1	3.1 ± 0.2	3.7 ± 0.2
		Yes	10.1 ± 0.8	4.4 ± 0.4	4.3 ± 0.3	4.3 ± 0.1	4.7 ± 0.3
	Dataset	Motion	No correction	Ratio	Generative	B1+ per contrast	Ratio & B1+ per contrast	Generative & B1+ per contrast
	#1	No	4.0 ± 0.3	4.1 ± 0.4	4.1 ± 0.3	4.1 ± 0.2	4.1 ± 0.4	4.1 ± 0.3
		Yes	5.8 ± 0.7	4.8 ± 0.3	4.9 ± 0.4	4.9 ± 0.3	3.8 ± 0.3	3.9 ± 0.3
7T	#2	No	4.8 ± 0.2	4.9 ± 0.3	4.9 ± 0.3	4.8 ± 0.2	5.0 ± 0.2	4.9 ± 0.2
		Yes	8.4 ± 0.4	6.1 ± 0.3	6.2 ± 0.3	8.0 ± 0.4	5.3 ± 0.1	5.5 ± 0.1
	#3	No	5.1 ± 0.2	5.0 ± 0.2	5.0 ± 0.2	5.0 ± 0.2	5.1 ± 0.3	5.1 ± 0.2
		Yes	13.4 ± 1.9	7.9 ± 1.1	8.3 ± 1.2	11.8 ± 0.9	6.3 ± 0.3	6.8 ± 0.3

FIGURE 5

Combining receive sensitivity and correction at 7T, for participant 3 (last row in Figure 4). The first row shows an example map without and with inter‐scan motion, before and after net receive sensitivity correction, and employing a separate map for each contrast in the case of inter‐scan motion. The second row shows (normalised) error maps with respect to the ground‐truth map . In this example, inter‐scan motion biases were corrected with the generative modeling approach. Each map combines a minimum of three scans (two VFA and one acquisitions). “ per contrast” correction incorporates one additional scan so as to map the transmit field at each position “Receive sensitivity correction” incorporates two (one per position) additional low‐resolution, single echo calibration scans (see section 2.2).

3T Validation

The net motion is summarised as the root‐sum‐of‐squares, across the three orthogonal axes, of the translations or rotations independently. The net translational and rotational motion in the “no motion” condition was 0.8 mm and 0.3°. These were increased to 1.2 mm and 18.1° in the inter‐scan motion case. When the T1w and PDw volumes were acquired in the same position, the MAE captured the test–retest variability, which was approximately 3% at 3T and 4%–5% at 7T. In the absence of overt motion, correcting for the differential sensitivity modulation did not substantially change the MAE. In the presence of overt motion, the MAE rose to 10%. It was reduced to 4.7% by the method of Papp et al. and to less than 4.4% by the proposed correction schemes, with or without incorporating the body coil in the generative modeling approach. The histograms in Figure 3 confirm that the method did not introduce any bias to the estimates.

Extension to 7T

At 7T, the range of motion varied across participants. The net translational and rotational motion in the “no motion” conditions did not exceed 1.6 mm and 1.0°, respectively. In the inter‐scan motion cases, the net translation ranged from 1.6 to 7.9 mm, while the net rotation ranged from 2.7 to 11.2°. Rotational motion led to more apparent artifacts. The overall amplitude of motion dictated the increase in MAE, which reached a maximum of 13.4% under the tested conditions (cf, motion summaries in Figure 4 and MAE in Table 1). The proposed correction scheme reduced the MAE (5%–8%), though not to the level of no overt motion.

The variability of the transmit field, , across head positions was found to be much higher than at 3T. Incorporating position‐specific maps reduced the MAE (5%–12%) even without correcting for the differential receive sensitivity modulation.

The greatest reductions in MAE were obtained by correcting for position‐specific transmit and receive fields, reaching 4%–7%, converging on the level obtained in the absence of overt motion (i.e. 4%–5%).

Comparison of methods

Overall, the ratio and generative modeling approaches to correcting the effects of differential relative sensitivities in maps performed similarly. The MAE was marginally lower for the generative modeling approach at 3T (0.1%) and for the ratio approach at 7T (0.5%). However, these differences were small relative to the variability across cases (Table 1).

Numerical error

At 3T the relative transmit efficiency ranged from 0.97 to 1.04, whereas the relative receive field (measured via ) ranged from 0.84 to 1.18. At 7T the relative transmit efficiency ranged from 0.85 to 1.18 under comparable motion conditions. These ranges were used in the simulations which revealed that without correction, inter‐scan motion caused error as high as 130%. Over the 4D parameter space investigated, a median of 29% of the error was caused by transmit field effects and 71% by receive field effects. Figure 6 shows a plane of this error as the relative transmit and receive fields change. Position‐specific offers only partial correction (Figure 6B). Larger error reduction arises from receive field correction (Figure 6C). Combining both (Figure 6D) shows receive field effects are removed but transmit sensitivity remains (when ).

FIGURE 6

error (in percentage of the true ) as a function of the relative transmit field and relative receive sensitivity between two head positions. Isocontours (black lines) and their associated percentage error are overlaid on each subplot. The four panels show this error with different degrees of correction: (A) none, (B) correction for position‐specific transmit field, (C) correction for position‐specific receive sensitivity, (D) both corrections. Note that a small amount of error remains, even with both corrections, because position‐specific transmit field effects lead to inaccuracies in the estimation of via κ 1,2, i.e.

DISCUSSION

We have introduced methods for correcting inter‐scan motion artfacts in quantitative MRI that do not rely on the availability or spatial homogeneity of a body coil. The approaches are based on estimating the relative sensitivity modulation across positions, and successfully reduced error in maps at both 3T and 7T.

At 3T, the proposed approaches outperformed a previously established correction method. This can be attributed to the fact that the method of Papp et al. assumes that the reference modulation of the body coil is independent of position, whereas the proposed methods do not. Instead they specifically account for the relative sensitivity across positions thereby restoring consistent modulations.

In the motion conditions tested here, both proposed approaches (ratio with Gaussian smoothing, or generative modeling) produced comparable improvements in reproducibility in the presence of inter‐scan motion. Equally importantly, when there was no overt motion neither method decreased reproducibility, which was at a level in keeping with previous reports for similar resolution MPM data. ,

The ratio method benefits from its simplicity, but may be vulnerable to low signal‐to‐noise ratio given that it defines one calibration image as the reference (denominator in equation (5)). The alternative generative modeling approach has the benefit of inherently adapting to variable signal‐to‐noise ratio by estimating the position‐specific net sensitivity modulation relative to a common image, which is their barycenter mean. This common image dictates the final modulation of all the corrected volumes. The generative model can also easily incorporate any additional data, for example, body coil images as done at 3T, which further flattens the final modulation. Furthermore, rigid registration could be interleaved with model fitting to reach a better global optimum. Finally, the generative model could naturally be integrated with any fitting approach that defines a joint probability over all acquired data, such as Balbastre et al. in the context of MPM.

The impact of movement on the effective transmit field has previously been investigated in the context of specific absorption rate management. , , , An important additional finding of the present work is the impact this can have on estimates at 7T, which was negligible at 3T as demonstrated previously.

Limitations

These methods are specifically designed for the correction of inter‐scan motion and therefore cannot address intrascan motion, which may be more likely to occur coincidentally with inter‐scan motion, for example, with uncompliant participants. Although the dominant source of error in was related to receive field effects, the MAE was further reduced by additionally accounting for the positional‐dependence of the transmit field. However, acquiring a map at multiple positions comes at a cost of increased scan time and inevitably leads to a greater temporal separation between the calibration data and those volumes it is used to correct. Issues such as this, coupled with other uncorrected effects, for example, position‐dependent effects (no reshimming was performed during the experiments), may underlie the fact that the corrections implemented do not reduce the MAE quite to the level of no motion. This finding recapitulates those of Papp et al., though the discrepancy is lower in this work, which is likely because the assumption of a flat body coil receive sensitivity is no longer made. The fact that even with combined receive and transmit field corrections, the MAE is never reduced to the level of no motion is in line with the simulations, which show that position‐dependent transmit field effects remain in the calibration data and propagate into the estimates.

An additional limitation of the generative model is its reliance on a Gaussian noise assumption, which is violated in the background (but not in the tissue, given the high signal‐to‐noise ratio of the calibration scans). Although we did not find this violation to hamper sensitivity estimation in the present study, the model could nonetheless be modified to incorporate a Rice or noncentral Chi likelihood.

CONCLUSIONS

Inter‐scan motion causes serially acquired weighted volumes to be differentially modulated by position‐specific coil sensitivities leading to substantial errors when they are combined to compute quantitative metrics. We have demonstrated the efficacy of two novel methods at reducing these artifacts in the context of mapping. The proposed methods do not require a body coil making them ideally suited for use at 7T, and can also be applied to the computation of other quantitative metrics, such as magnetization transfer saturation, that similarly assume constant modulation across multiple weighted acquisitions. Given that the acquisition of the receive field calibration data is rapid, and their use for receive field correction does not degrade reproducibility in the no inter‐scan motion condition, we would recommend that this correction routinely be incorporated into qMRI workflows. Although specifically at 7T motion will always induce a degree of transmit field change and associated error, accurate transmit field calibration is prohibitively long for routine repetition. Future work will therefore explore more time efficient approaches of correcting this residual error.

CONFLICT OF INTEREST

The Wellcome Centre for Human Neuroimaging have institutional research agreements with Siemens Healthcare.