Time-resolved angiography with stochastic trajectories for dynamic contrast-enhanced MRI in head and neck cancer: Are pharmacokinetic parameters affected?

Purpose: To investigate the e ff ects of di ff erent time-resolved angiography with stochastic trajectories (TWIST) k -space undersampling schemes on calculated pharmacokinetic dynamic contrast-enhanced (DCE) vascular parameters. Methods: A digital perfusion phantom was employed to simulate e ff ects of TWIST on characteristics of signal changes in DCE. Furthermore, DCE-MRI was acquired without undersampling in a group of patients with head and neck squamous cell carcinoma and used to simulate a range of TWIST schemes. Errors were calculated as di ff erences between reference and TWIST-simulated DCE parameters. Parametrical error maps were used to display the averaged results from all tumors. Results: For a relatively wide range of undersampling schemes, errors in pharmacokinetic parameters due to TWIST were under 10% for the volume transfer constant, K trans , and total extracellular extravascular space volume, V e . TWIST induced errors in the total blood plasma volume, V p , were the largest observed, and these were inversely dependent on the area of the fully sampled k -space. The magnitudes of errors were not correlated with K trans , V p and weakly correlated with V e . Conclusions: The authors demonstrated methods to validate and optimize k -space view-sharing techniques for pharmacokinetic DCE studies using a range of clinically relevant spatial and temporal patient derived data. The authors found a range of undersampling patterns for which the TWIST sequence can be reliably used in pharmacokinetic DCE-MRI. The parameter maps created in the study can help to make a decision between temporal and spatial resolution demands and the quality of enhancement curve characterization. C 2016 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License.


INTRODUCTION
View-sharing techniques allow for significant reduction of acquisition time, enabling fast high-resolution measurements without compromising anatomical coverage. Such techniques are often used in cardiac function examinations, 1,2 angiography, 3 flow and perfusion imaging, [4][5][6][7] and situations where breath-holds are required. 8,9 View-sharing methods employ partial k-space updating with assumptions of a priori knowledge to recover spatial resolution [10][11][12][13] or coverage loss. 14 Usually, it is possible to vary the density of sampling in order to match acquisition parameters with clinical temporal and spatial redundancy (e.g., respiratory or cardiac motion and vascular input function), and thus, reduce errors. In the case of time-resolved angiography with stochastic trajectories (TWIST), 11 it is possible to vary the size of the central (fully sampled) and peripheral parts of k-space, together with the percentage of periphery sampled at a given acquisition time-point. In the case of 4D acquisition, the use of a small central compartment results in increased temporal resolution. This is especially beneficial for dynamic contrastenhanced (DCE) MRI, where a high temporal resolution is required to appropriately characterize a progressive signal intensity change after injecting gadolinium-based (Gd-based) intravenous contrast. 15 Typical temporal resolution in these types of measurements is between 2 and 10 s depending on the level of perfusion, kinetic model used, and the necessity of motion compensation (i.e., breath holding or cardiac triggering). Variable flip angle 4D gradient echo sequences are commonly used in DCE-MRI to calculate longitudinal relaxation time, T 1 , using a signal ratio between the proton density (low flip angle) and T 1 -weighted (higher flip angle) images. T 1 relaxation time changes allow for Gd concentration calculations as a function of time. 16 The series of T 1 -weighted images is usually acquired for a period of 3-10 min, allowing the characterization of uptake and washout of the contrast agent in studied tissues. Using appropriate kinetic models, the dynamics of Gd concentration can be related to vascular parameters such as the volume transfer constant between blood plasma and extracellular extravascular space (K trans ), the total extracellular extravascular space volume (V e ), or the total blood plasma volume (V p ). View-sharing techniques can be compatible with these methods. However, undersampling of peripheral parts of k-space may lead to impaired characterization of enhancement curves, ringing artifacts, and impaired fat suppression. [17][18][19] In this study, we investigate the effects of different kspace undersampling schemes with TWIST on calculated pharmacokinetic (PK) DCE vascular parameters. High temporal resolution DCE data from a group of patients with head and neck squamous cell carcinoma (HNSCC) were used to evaluate TWIST-induced errors in K trans , V e , and V p .

METHODS
MR images acquired prior to treatment in eight patients with HNSCC were used in the study. Written informed consent was obtained from all patients in this study, which was approved by the institutional review board (Royal Marsden Hospital Committee for Clinical Research) and research ethics committee (NHS REC number 10/H0801/32).

2.A. MRI protocol
DCE data acquired with a high temporal resolution without the use of TWIST were used as the ground truth reference. MRI was performed at 1.5 T (Philips Intera, Philips Medical Systems, Best, Netherlands). Patients were imaged in a 5-point thermoplastic mask and two flexible surface coils centered over the volume of interest were used. The DCE protocol included a trans-axial 3D spoiled fast gradient echo sequence (TE = 1 ms, TR = 4 ms, FOV = 256 × 256 mm 2 , acquisition/reconstruction matrix: 112 × 128/128 × 128, 10 slices, with 2 × 2 × 5 mm 3 voxels, parallel imaging factor of 1.7, and partial Fourier acquisition of 60% in the anterior-posterior direction). A series of 20 proton density-weighted images (flip angle, FA = 4 • ) was initially acquired prior to contrast injection. This was followed by 100 T 1 -weighted acquisitions (FA = 20 • ) obtained sequentially with 1.5 s temporal resolution. Gd contrast was injected intravenously at the start of the tenth dynamic scan as bolus through a peripherally placed cannula using an automatic injector (0.2 ml/kg body mass, 3 ml/s injection rate, Dotarem, Guerbet, France) and followed by a saline flush (20 ml at 3 ml/s).

2.B. PK modeling
Regions of interest (ROIs) were manually delineated by a radiation oncologist together with an expert MRI radiologist on T 1 -weighted images at the time-point with maximum contrast concentration (e.g., 11th dynamic measurement, approximately 15 s after the start of contrast injection) for each slice, around all primary tumors and involved lymph nodes. Signal ratio between the proton density (FA = 4 • ) and T 1 -weighted (FA = 20 • ) images was used to calculate T 1 relaxation time. First five images in each dynamic acquisition were excluded from analysis to ensure magnetization equilibrium. Signal from remaining 15 proton-weighted images acquired before the contrast injection was averaged in order to reduce the potential influence of incidental motion (i.e., due to swallowing). Changes of T 1 were used for Gd concentration calculations as a function of time, used in the PK modeling. The gadolinium contrast onset (∼10th dynamic) time was manually adjusted for each DCE session. The data were analyzed using the software package MRIW (Institute of Cancer Research, UK) 20 with the extended Kety model 21 and a population-based arterial input function (AIF). 22 A set of parameters was derived including: K trans , V e , and V p . DCE maps were produced for each parameter and lesion ROI.

2.C. Numerical DCE phantom
A digital perfusion phantom 23 was employed to simulate effects of TWIST undersampling on characteristics of signal changes in DCE-MRI using  (Mathworks, Cambridge MA). The extended Kety model was used to generate enhancement curves for K trans in a range of 0.125-0.7 min −1 (fixed V e = 0.25, V p = 0.01), V e in a range of 0.125-0.7 (K trans = 0.25, V p = 0.01), V p in a range of 0.03-0.2 (K trans = 0.25, V e = 0.25), and temporal resolution of 0.5 s. Simulated homogenous lesions (radius = 3 pixels) were positioned along the phantom perimeter and a set of 300 simulated images (128 × 128 pixels) was generated.
F. 1. TWIST acquisition scheme for A = 50% and B = 33% (P-phase encoding, F-frequency encoding, and S-partition/slice direction): k-space was divided into central (A 1 ) and 3 peripheral parts (B 1-3 ). In the first acquisition (t 1 ), the full k-space [A 1 (t 1 ), B 1 (t 1 ), B 2 (t 1 ), B 3 (t 1 )] was sampled. In the next acquisition, only the central A 1 (t 2 ) and 33% of peripheral B 1 (t 2 ) parts of the k-space were sampled and the missing data of B 2-3 (t 2 ) copied from adjacent time point t 1 : . Similarly, in the following acquisition steps, all points in compartment A and a percentage of the points in compartment B were sampled, with missing portions of B copied from adjacent time points.

2.D. TWIST simulations
A set of images measured or generated at each DCE timepoint was Fourier transformed to obtain k-space data using . Each k-space data set was divided into central and peripheral parts (Fig. 1). Five different areas of the central part (A = 2%, 10%, 20%, 33%, and 50%) were considered. The peripheral part was subsequently undersampled, such that a defined percentage of remaining k-space points (B = 10%, 33%, and 50%) was randomly sampled. For each image, all points in compartment A and a percentage of the points in compartment B were used, with missing portions of B copied from adjacent time points.
Created k-spaces were Fourier transformed to obtain a set of simulated TWIST images used for DCE pharmacokinetic calculations (Fig. 2). This approach is similar to that used F. 2. Schematic diagram of TWIST simulation: fully sampled magnitude image sets were Fourier transformed to obtain k-space data. Five different sizes of the central part (A = 2%, 10%, 20%, 33%, and 50%) were considered. The compartment A was fully sampled while B was sampled with reduced density (B = 10%, 33%, and 50%) with missing portions of B copied from preceding time points. Simulated k-spaces were Fourier transformed to obtain a set of TWIST images used for DCE calculations. by Song et al. in a simulated renography phantom study. 17 Ratios of TWIST undersampled and fully sampled reference acquisition times (TA TW ) were calculated for different TWIST combinations.

2.E. Error calculation for DCE parameters
In order to determine the magnitude of K trans , V e , and V p errors, sets of DCE TWIST images were simulated for the digital phantom and every patient, with various TWIST parameter combinations. Absolute errors were calculated as a difference between measured fully sampled reference (DCE ref ) and TWIST-simulated (DCE TWIST ) DCE parameters on a pixel-by-pixel basis, where x and y are given pixel coordinates. The percentage errors were than calculated as follows: The median percentage difference was calculated for all voxels within each patient ROI (n = 8) and for 15 simulated A/B combinations (A = 2%, 10%, 20%, 33%, 50% and B = 10%, 33%, 50%). In the case of the phantom simulations, 12 TWIST combinations (A = 10%, 20%, 33%, 50% and B = 20%, 33%, 50%) were used and the results were presented in Table I. In the next step, the TWIST error distribution was studied for the whole group of patients. The results obtained for patient ROIs were averaged for each undersampling pattern and the results also presented in the form of parametrical TWIST error maps. Standard deviations (SDs) of the TWIST parameter errors were calculated and also reported in the form of parametrical maps, allowing for the assessment of intrapatient variability of parameter errors.
Kendall's tau (τ) 24 was used to test for correlations between K trans , V e , V p , and corresponding percentage errors for all voxels. The null hypothesis stated that there was no correlation between measured DCE parameters and percentage errors.
The strength of correlation was tested 24 and was considered significant if p < 0.01. Statistical analysis was performed using the  Statistics Toolbox.

3.A. Phantom simulations
The effect of TWIST on DCE parametrical maps simulated for A,B = 20%,20% is presented in Fig. 3. Relatively high kspace undersampling results in blurring and ringing artifacts. DCE TWIST error maps are presented in Table I. The error magnitudes depended on the choice of undersampling pattern, parameter type, and its value. The highest errors were observed for V p (max error 30%) and the smallest for K trans , with errors below 15% for all measured values. There was a pronounced increase of error for A < 33%, whereby the error increased with the value of the DCE parameter.

3.B. Clinical simulations
Patient characteristics and tumor stages are summarized in Table II. The results for TWIST DCE parameter errors are  presented in Table III Median TWIST percentage error maps (HNSCC patients) are presented in Fig. 4. The K trans and V e differences were under 10% for a relatively wide range of undersampling schemes: A > 30% and B > 20%. TWIST had the largest effect on the V p parameter, where differences lower than 15% were only achieved using less extensive view-sharing parameters of A > 35% with an extended region between B: 25-35 and A > 25%.
Intrapatient variability of the K trans error was greater for small A and B with error standard deviation >5% for TWIST parameters A and B < 25%. Standard deviations of total extracellular extravascular space volume errors were generally higher (max SD = 10%) than K trans for all TWIST schemes. The variability of plasma volume errors was higher for small B values with a steep increase at B < 20%. There was no significant correlation (median τ = −0.12, p = 0.09) observed between the K trans values calculated for all measured voxels and corresponding TWIST percentage errors for all undersampling regimes used. Similarly, there was no correlation for V p (median τ = 0.09, p = 0.11) and a weak, statistically significant, negative correlation (median τ = −0.16, p < 0.001) for V e parameters and measured percentage errors. There was also no correlation between tumor ROI size and TWIST percentage errors (K trans τ = 0.35, p = 0.28; V e : τ = −0.29, p = 0.4, V p : τ = 0.14, p = 0.55).

DISCUSSION
The high temporal resolution of the conventional 3D gradient echo sequence allows for reliable calculations of DCE parameters. However, short acquisition times are achieved at a cost of decreased spatial resolution or volume coverage, which can lead to significant partial volume effects. This can be problematic in the case of longitudinal DCE measurements during a course of treatment, where decreasing tumor volume in responding lesions may cause the number of viable tumor voxels to become limited. Limited superior/inferior lesion coverage can also be troublesome in HNSCC where local T III. TWIST DCE parameter errors, ROI sizes, and acquisition time ratios (TA TW ) for all patients and undersampling combinations.
A,B (%) 50,50 50,33 50, 10  lymph nodes are a common site of metastatic disease and should be imaged together with the primary site. In the case of HNSCC, the desirable superior/inferior coverage is between 10 and 20 cm with a 20-25 cm in-plane field of view, allowing for full head and neck anatomy coverage. There is increasing interest in the use of functional imaging for biological target volume delineation in radiotherapy treatment planning. 25,26 Identification of less well-perfused regions of the tumor is of interest, as these may reflect hypoxic subvolumes, and therefore, areas of relative radioresistance. 27,28 However, this requires not only in-plane but also through-plane high spatial resolution in order to accurately identify and delineate biological tumor subvolumes, which could be potentially used for a local dose boosting. Confidence in the coregistration of radiotherapy dosimetric data and DCE parameters is particularly important where steep dose gradients exist between target volumes and surrounding healthy tissue. This could be further improved, thanks to the recent development of MR-based RT planning 29 and development of MR-Linac. 30 The use of TWIST can overcome these limitations of conventional 3D gradient-echo techniques by increasing the user's ability to find a good compromise between coverage, temporal, and spatial resolution.
We demonstrated the influence of the TWIST view-sharing technique on DCE calculations using a digital phantom and heterogeneous groups of primary and nodal tumor sites. The use of simulated phantom data enabled the generation of reference DCE data not influenced by a limited temporal resolution and for a wide range of vascular parameters (K trans , V e , and V p ). We found that the plasma volume was the most significantly affected parameter by the TWIST undersampling, which can be explained by a dependence of V p on the initial Gd uptake peak. K trans and V e parameter errors were below 10% for a wide range of TWIST combinations and K trans and V e values, as presented in Table I. The data F. 4. TWIST parameter error maps-clinical results. An absolute value of mean percentage difference of K trans , V e , and V p calculated using reference (no TWIST) and simulated TWIST DCE image sets is shown on the left. The white iso-contours indicate areas of all K trans , V e , and V p errors below 15%. Standard deviations for the percentage differences are shown on the right. Simulated TWIST undersampling partition sizes: A = 2%, 10%, 20%, 33%, 50% and B = 10%, 33%, 50%).
suggest that it is possible to achieve errors below 15% in all K trans , V e , and V p parameters for A > 20% and V p below 0.1, which is expected in HNSCC (median V p = 0.022 ± 0.012 in this study).
Error maps created for various undersampling schemes can help to assist in the choice of TWIST parameters to be applied for DCE studies. For example, it would be possible to improve the spatial resolution of the DCE experimental protocol used in this study from 2×2×5 mm 3 to 2×2×2.75 mm 3 employing A,B = 33%,33% TWIST undersampling, maintaining high temporal resolution (1.5 s) and coverage. In this case, the predicted TWIST induced errors would be below 10% for K trans , V e and under 15% for the V p . The resulting spatial resolution is in line with routine radiotherapy treatment plans taking into consideration patient positioning, motion-related uncertainties, and deliverable treatment beam geometry. 31 The potential detrimental effects of view-sharing methods on MRI image quality were recognized before. Most studies however focused on qualitative assessment of clinical images [17][18][19] with few performing a quantitative analysis of dynamic signal errors. 32,33 The majority of these studies employed the "keyhole" method, which is an extreme case of view-sharing, where only the centre of k-space is updated continually, while the periphery of k-space is acquired only once at the beginning, and optionally again at the end of the dynamic series. 12 This leads to a strong geometrical dependence of dynamic signal errors and the requirement for the central kspace size to be restricted by the approximate minimum size of the expected lesion. 32 The TWIST method is less sensitive to these effects due to a sparse but continuous acquisition of peripheral part of the k-space. In our study, we do not observe correlation between lesion size and DCE TWIST errors.
One of the requirements of DCE-MRI is a requirement for a reliable measure of the AIF, representing a time course of the contrast concentration in the blood plasma pool. Patient-specific measurement of the AIF from a vessel feeding the tumor would be desirable. In practice, however, a population-based AIF is commonly used for measuring cohort and individual vascular parameters since the patient-specific repeatability is hampered by in-flow effects, nonlinear signal response, nonuniform B 1 field, pulsatile flow, and partial volume effects. 34 This approach was also adapted in this study. However, it is expected that higher temporal resolution achieved with TWIST could be beneficial in characterization of patient-specific AIF. Song et al. 17 found A,B = 20%,20% TWIST sampling optimal in characterization of simulated normal and diseased kidney function (glomerular filtration rate and renal plasma flow errors <10%) with one aortic input function.
Our work has several limitations. First, the k-space is obtained by Fourier transforming magnitude images, and as a consequence, phase effects on the signal are disregarded. This simulation applies either to situations where field inhomogeneity is negligible or to situations where the central portion of the k-space is sufficiently sampled to map field inhomogeneity. In a real TWIST-DCE experiment, the 3D k-space is sparsely sampled in the phase and compartment (slice) directions (Fig. 1). In our work, data are effectively undersampled in phase and readout direction and therefore is an equivalent of sagittal or coronal rather than axial 3D acquisition. Another limitation is the fact that the sampling frequency of the simulated data was constant and determined by the acquired reference DCE data. This therefore represents the use of TWIST to increase spatial resolution or coverage. This is not the case if the TWIST is used to improve temporal resolution where the sampling frequency would be higher and the induced DCE errors are likely to be lower than presented in our study. Finally, our limited patient cohort size does not represent full range of possible vascular and spatial tumor characteristics and our clinical evaluation is specific to HNSCC tumors. This was addressed to a certain extent using the digital phantom, which enabled to study TWIST effects for a wide range of DCE parameters and without temporal resolution limitations. In our work, we combined a DCE phantom with the well-known geometry (Shepp-Logan) 23 and Qualitative Imaging Biomarkers Alliance (QIBA). 35 However, a disadvantage of this approach is an arbitrary geometry of simulated enhancing lesions, the lack of spatial heterogeneity, and fixed size of regions of interest. As a consequence, the TWIST error magnitude can be elevated for particular undersampling patterns, for which effects of spatial blurring and compromised edge encoding populate significant subvolumes of simulated ROIs. This could explain higher errors calculated with A = 20%, when compared to A = 10%. The main advantage of using clinical data from Head & Neck studies is to reproduce accurately the range of spatial frequencies present in the region of interest, and thus evaluate more accurately the error introduced into pharmacokinetic parameters by the use of view sharing techniques. It also allows for an appropriate level of physiological noise influencing the quality of pharmacokinetic modeling. Despite of these shortcomings, our simulation achieves its objectives, as it eliminates an appropriate amount of information in space and time leading to distortion caused by TWIST to the PK parameters in HNSCC mapped. In our clinical simulations, we did not observe correlations between K trans , V e , V p values, and the magnitude of error, for a heterogeneous group of tumor ROIs. This is in line with the phantom simulations showing pronounced increase of error for DCE parameters higher than those observed in the studied HNSCC patient group. Furthermore our simulations provide an evaluation of the minimum sampling requirements to keep the error level on the determination of pharmacokinetic parameters K trans , V e , and V p under a given threshold. In the case of HNSCC tumors (median K trans = 0.194 min −1 , V e = 0.223, and V p = 0.022 in this study), the recommended size of the TWIST undersampling is A and B > 20%. This allows achieving errors below 15% in all DCE parameters.

CONCLUSION
In conclusion, we demonstrated a method to validate and optimize k-space view-sharing techniques used for pharmacokinetic DCE studies in head and neck cancers. In this setting, the TWIST sequence can be used reliably for pharmacokinetic DCE-MRI using a range of under-sampling patterns. The parameter maps created in the study can be used to balance temporal and spatial resolution demands to allow optimal enhancement curve characterization.