Estimation of contrast agent bolus arrival delays for improved reproducibility of liver DCE MRI

Local ethics committee approval was obtained and participants provided informed written consent. Volunteers were recruited via advertisement within the university campus and were eligible if (a) they had no MRI contraindication, (b) were not taking any long-term medication (excluding the oral contraceptive pill) and (c) had no documented history of previous liver or gastrointestinal disease. Fourteen volunteers were screened of which one was excluded because of claustrophobia. The final cohort consisted of seven males (aged 26.5  ±  1.4 years) and six females (aged 31.2  ±  2.6 years). Participants fasted for 6 h prior to MRI and avoided caffeinated fluids. A 19G cannula was sited in a peripheral upper limb vein in preparation for administration of contrast. The breathing protocol was then explained to subjects before entering the scanner by the study coordinator (radiology research fellow with 5 years experience). For reproducibility studies, nine subjects consented to be re-scanned 7 days after the original study following identical preparation and MRI protocol, at a comparable time of the day (within 2 h).

Imaging was performed using a 3.0T scanner (Achieva, Philips Healthcare, Best, Netherlands) using a 16 channel body coil (SENSE XL-Torso, Philips Healthcare, Best, Netherlands). After initial anatomical imaging using a breath hold balanced steady-state free precession (SSFP) sequence, DCE studies were planned to ensure inclusion of the whole liver volume, retroperitoneal great vessels and the heart. A multi-flip angle T1 measurement was undertaken using three-dimensional (3D) gradient echo imaging at five different flip angles (5, 7, 10, 15 and 20°), with phase based B₁ mapping for B₁ non-uniformity correction (Treier et al 2007). DCE imaging in the coronal plane to minimise inflow effects, was performed using a 3D gradient turbo field echo (TFE) imaging with spectral attenuated inversion recovery (SPAIR) fat suppression. Thirty overcontiguous slices were acquired, then interpolated to sixty, with a total dynamic scan time of 3.35 s per 15 cm volume, scanned sequentially for 5 min (sequence parameters given in table 1). After the first five volumes were acquired, ten ml of Gd-DOTA (gadoterate dimeglumine, Dotarem^®, Guerbet, Roissy, France), diluted in 10 ml of normal saline, was injected at 4 ml s⁻¹ (Spectris^®, Medrad Inc., USA), followed by a 20 ml saline flush. Subjects were given the first breath hold instruction before the CA injection: this minimised the likelihood of motion artefact in the early part of the DCE study (including the VIF peaks). Thereafter they were asked to continue self-directed breath holds in expiration for the duration of the study.

Post-processing was performed using Matlab code (MathWorks, Natick, USA) developed in house. DCE volumes corrupted by significant motion artefact noise were discarded (average of 21/90 volumes discarded). No VIF peaks were missed in discarded data sets. Five coronal slices centred around the portal vein, each separated by 10 mm were then selected for analysis. Each slice was matched to data from the previously derived T1 maps and registered using robust data decomposition registration to correct for tissue displacement and deformation (Hamy et al 2014). Missing SI data from discarded volumes was estimated using linear interpolation. The interval for interpolation only exceeded one discarded volume on two occasions, both well after VIF and parenchymal peak SIs. In these instances, interpolation for two successive discarded volumes took place. Pixel wise conversion of sequential post-contrast SI into CA concentration was then undertaken using previously described methods for each of the five slices (Aronhime et al 2014, Gill et al 2014). Three parenchymal ROIs were positioned on each slice (total 15 ROIs), firstly in the right upper region (segments VII/VIII), left liver (segments II/III) and right lower region (segments V/VI). Care was taken to ensure parenchymal ROIs excluded any major inflow or outflow vessels (HA, PV and hepatic venous radicles). ROIs were then also positioned within the left ventricle and PV to derive each VIF (figure 1), as the left ventricle in our experience delivered more consistent VIFs. Perfusion parameters (detailed below) extracted from all fifteen ROIs (three ROIs on five slices) were averaged across all subjects for different post-processing method comparisons. All post-processing was undertaken by the study coordinator (radiology research fellow with 5 years' experience of abdominal MRI).

Zoom In Zoom Out Reset image size

Dual input single compartment modelling was undertaken as reported previously (Materne et al 2002, Hagiwara et al 2008). Briefly, liver parenchymal CA concentration as a function of time (${{C}_{\text{L}}}(t)$) can be expressed as:

$$\begin{eqnarray}&&{{C}_{\text{L}}}(t)={\int}_{0}^{t}\left[{{k}_{1\text{a}}}{{C}_{\text{a}}}\left({{t}^{\prime}}-{{\tau}_{\text{a}}}\right)+{{k}_{1\text{p}}}{{C}_{\text{p}}}\left({{t}^{\prime}}-{{\tau}_{\text{p}}}\right)\right]{{\text{e}}^{-{{k}_{2}}\left(t-{{t}^{\prime}}\right)}} ~\text{d}{{t}^{\prime}}\,\end{eqnarray} \tag{ 1 }$$

where ${{C}_{\text{a}}}(t)$ represents the arterial input CA concentration as a function of time, ${{C}_{\text{p}}}(t)$ represents the PV input CA concentration as a function of time, ${{k}_{1\text{a}}}$ represents the arterial inflow constant, ${{k}_{1\text{p}}}$ represents the PV inflow constant, ${{k}_{2}}$ represents the outflow constant, ${{\tau}_{\text{a}}}$ represents the delay between the arrival of CA in the AIF and parenchymal ROIs and ${{\tau}_{\text{p}}}$ represents the delay between arrival of CA in the PVIF and parenchymal ROIs (figure 2(a)). Model fitting was undertaken using non-linear least squares fitting with in house developed Matlab code. Inflow and outflow constants were used to derive estimates of PV perfusion (ml/min/100 g), (TLBF, sum of HA and PV perfusion, ml/min/100 g), HA fraction (%), distribution volume (DV, %) and mean transit time (MTT, seconds) as reported previously (Materne et al 2002, Hagiwara et al 2008).

Zoom In Zoom Out Reset image size

CA bolus arrival delays in the modelling process effectively shift the VIFs in the modelling process forward by increments defined by the temporal resolution of the data (3.35 s in this study) (figure 2(a)). To investigate the effects of altering AIF and PVIF CA bolus arrival delays, pharmacokinetic modelling was undertaken after introducing successive increments to ${{\tau}_{\text{a}}}$ (6 steps, up to 20.10 s) and ${{\tau}_{\text{p}}}$ (3 steps, up to 10.05 s) for each dataset (larger delays were not studied as these are less likely to be physiological). Perfusion parameters for each combination of VIF delays (n  =  18) were then averaged across subjects for analysis.

Perfusion parameters and their 7 days reproducibility were then compared using the following methods for handling CA bolus arrival delays:

• (i)  
  No delaysAssuming no delay between VIFs and parenchymal enhancement (i.e. ${{\tau}_{\text{a}}}$ and ${{\tau}_{p}}$ both set to zero), fixed across all datasets (Murase et al 2007, Miyazaki et al 2008).
• (ii)  
  Freely modelled delaysFree (unconstrained) modelling of AIF and PVIF delays to optimise model fitting by minimising the residual sum of squares (Hagiwara et al 2008, Kim et al 2008).
• (iii)  
  Pre-estimated delays with constrained free modellingThe first five data points for VIF and parenchymal CA concentrations were used to determine the baseline pre-CA concentration for each curve. Linear regression between the first VIF data point exceeding the upper limit of the 95% confidence interval of baseline pre-CA concentration, and the VIF peak was undertaken for the AIF and PVIF. The VIF CA bolus arrival times (${{t}_{{{C}_{\text{a}}}\text{arrival}}}$ and ${{t}_{{{C}_{\text{p}}}\text{arrival}}}$) were then each estimated as the point of intercept of the regression line with the time axis. Parenchymal enhancement was less noisy because of larger ROI size and lower susceptibility to flow artefact. The parenchymal CA arrival time (${{t}_{{{C}_{\text{L}}}\text{arrival}}}$) was then defined as the last data point before the parenchymal CA concentration exceeded the upper limit of the 95% confidence interval of the baseline pre-CA concentration. Figure 2(b) demonstrates the process in detail. Pre-estimates for ${{\tau}_{\text{a}}}$ and ${{\tau}_{\text{p}}}$ were then determined as:

  $$\begin{eqnarray}&&\tau _{\text{a}}^{\prime}={{t}_{{{C}_{\text{L}}}\text{arrival}}}-{{t}_{{{C}_{\text{a}}}\text{arrival}}}\end{eqnarray} \tag{ 2 }$$

  $$\begin{eqnarray}&&\tau _{\text{p}}^{\prime}={{t}_{{{C}_{\text{L}}}\text{arrival}}}-{{t}_{{{C}_{\text{p}}}\text{arrival}}}\end{eqnarray} \tag{ 3 }$$

  As $\tau _{\text{a}}^{\prime}$ and $\tau _{\text{p}}^{\prime}$ represented estimates of VIF delays, limited by temporal resolution (3.35 s), the pre-estimates were then used to constrain the range in which pharmacokinetic free modelling of ${{\tau}_{\text{a}}}$ and ${{\tau}_{\text{p}}}$ could occur, to one time point before and one time point after each estimate (i.e. within a 6.7 s window).

To investigate the effect of altering AIF and PVIF CA bolus delays, calculated perfusion parameters for each delay were expressed as a percentage of those obtained when assuming zero delay between VIFs and parenchymal enhancement.

Kolmogorov–Smirnov tests were then used to confirm the normality of perfusion parameters derived (i) assuming no VIF delays (i.e. ${{\tau}_{\text{a}}}$ and ${{\tau}_{\text{p}}}$ both set to zero), (ii) freely modelled delays and (iii) pre-estimated delays with constrained free modelling. Repeated measures one-way analysis of variance (ANOVA) with corrections for non-sphericity were used to compare perfusion parameters using each of the three approaches to VIF delay estimation. Post hoc Tukey's test was then applied where significant differences were identified. Where variables were found not to be normally distributed, the Kruskal–Wallis test was used followed by post hoc Dunn's test if significant differences were identified. Paired t-tests/Wilcoxon matched-pairs signed rank tests as appropriate were used to compare VIF delay estimations obtained using freely modelled and pre-estimated delays with constrained free modelling. Seven-day reproducibility (n  =  9) was assessed using Bland–Altman (BA) analysis of agreement, with calculation of the mean difference (bias), 95% limits of agreement (LoA) and coefficients of variation. The threshold of statistical significance was defined to be p  <  0.05.

Calculated perfusion parameters for each AIF and PVIF delay increment expressed as a percentage of those obtained when assuming zero delay (i.e. ${{\tau}_{\text{a}}}$ and ${{\tau}_{\text{p}}}$ set to zero) are shown in figure 3.

Zoom In Zoom Out Reset image size

Estimated PV perfusion (figure 3(a)) decreased by as much as 31% (13.40 s AIF CA bolus arrival delay). Introducing PVIF CA bolus arrival delays increased PV perfusion by as much as 30% (10.05 s delay). A similar trend was demonstrated for TLBF (figure 3(b)), with perfusion estimates decreasing by as much as 10% (10.05 s AIF CA bolus arrival delay). Introducing PVIF CA bolus arrival delays increased TLBF by as much as 43% (10.05 s delay).

Because of small HA fraction estimates obtained when assuming zero VIF CA bolus arrival delays, introduction of CA bolus delays for both AIF and PVIFs resulted in increases of as much as 3247% (13.40 s AIF delay, 10.05 s PVIF delay, figure 3(c)). MTT increased by as much as 150% (0 s AIF delay, 10.05 s PVIF delay, figure 3(d)). DV reduced by up to 10% (0 s and 20.10 s AIF delay, 10.05 s PVIF delay, figure 3(e)).

Perfusion parameters calculated using each of the three methods are shown in table 2 and graphically in figure 4. HA fraction, DV and 'τ_p' were not normally distributed and underwent non-parametric statistical testing. Significant differences were demonstrated between the three methods for PV perfusion (F(1.25, 24.96)  =  7.29; p  =  0.0085) and HA fraction (H  =  23.94; p  <  0.0001), but there were no significant differences for the other parameters. Post hoc testing demonstrated significant differences between zero CA bolus VIF arrival delays and freely modelled delays for both PV perfusion and HA fraction, but only for HA fraction when free modelling was constrained using pre-estimated delays. Significant differences were demonstrated between freely modelled and pre-estimated delays with constrained free modelling for AIF (mean difference  −3.7  ±  1.1 s, p  =  0.0035) and PVIF CA bolus arrival delays (median difference 1.12 s, p  =  0.029).

Zoom In Zoom Out Reset image size

Reproducibility was assessed in 9 normal volunteers 7 days after the initial study (table 3, figures 5 and 6). The mean difference and BA 95% LoAs for repeated PV perfusion and TLBF measurements were smallest using pre-estimated delays with constrained free modelling. The coefficient of variation using this method was similar to unconstrained free modelling for both parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mean difference and BA 95% LoAs between repeated HA fraction measurements was smallest using no delays. The coefficient of variation was smallest using freely modelled delays. The smallest mean difference for repeated MTT and DV measurements was demonstrated using freely modelled delays. BA 95% LoAs and coefficients of variation were similar across all three methods for both MTT and DV (table 3).