Volumetric coronary endothelial function assessment: a feasibility study exploiting stack‐of‐stars 3D cine MRI and image‐based respiratory self‐gating

Abnormal coronary endothelial function (CEF), manifesting as depressed vasoreactive responses to endothelial‐specific stressors, occurs early in atherosclerosis, independently predicts cardiovascular events, and responds to cardioprotective interventions. CEF is spatially heterogeneous along a coronary artery in patients with atherosclerosis, and thus recently developed and tested non‐invasive 2D MRI techniques to measure CEF may not capture the extent of changes in CEF in a given coronary artery. The purpose of this study was to develop and test the first volumetric coronary 3D MRI cine method for assessing CEF along the proximal and mid‐coronary arteries with isotropic spatial resolution and in free‐breathing. This approach, called 3D‐Stars, combines a 6 min continuous, untriggered golden‐angle stack‐of‐stars acquisition with a novel image‐based respiratory self‐gating method and cardiac and respiratory motion‐resolved reconstruction. The proposed respiratory self‐gating method agreed well with respiratory bellows and center‐of‐k‐space methods. In healthy subjects, 3D‐Stars vessel sharpness was non‐significantly different from that by conventional 2D radial in proximal segments, albeit lower in mid‐portions. Importantly, 3D‐Stars detected normal vasodilatation of the right coronary artery in response to endothelial‐dependent isometric handgrip stress in healthy subjects. Coronary artery cross‐sectional areas measured using 3D‐Stars were similar to those from 2D radial MRI when similar thresholding was used. In conclusion, 3D‐Stars offers good image quality and shows feasibility for non‐invasively studying vasoreactivity‐related lumen area changes along the proximal coronary artery in 3D during free‐breathing.


| INTRODUCTION
Coronary atherosclerosis and its progression are spatially heterogeneous processes responsible for most cardiovascular events. 1 While healthy coronary artery endothelial cells in response to certain stressors 2 release protective substances such as nitric oxide (NO) that are anti-thrombotic, antiinflammatory and increase blood flow and arterial diameter, 3 abnormal endothelial cells are characterized by reduced NO release, and attenuated blood flow and arterial dilatation responses. This latter phenomenon of abnormal coronary endothelial function (CEF) is important because it is a marker of early local atherosclerosis, a predictor of sites of future atherosclerotic progression, an independent predictor of cardiovascular events and a responder to interventions shown to reduce cardiovascular risk. [3][4][5][6][7][8] Thus, CEF is often considered an important "barometer" of vascular health. 9 Conventional angiographic methods for CEF assessment typically use local administration of intracoronary endothelial-dependent stressors, cold pressor testing or isometric handgrip exercise (IHE) 3,10 to elicit impaired vasoreactivity, which predicts subsequent coronary artery disease (CAD) progression. 4 While catheterization-based methods enable assessment of the spatial heterogeneity of vasoreactivity response along a coronary artery, 3 their invasive nature precludes CEF assessment in healthy, low-risk and stable populations and thus has limited the impact and use of CEF measures in clinical practice and research. Non-invasive MRI measures of CEF 11 were introduced with electrocardiogram (ECG)-triggered two-dimensional (2D) spiral cine MRI and IHE, 12 an endothelial-dependent stressor. 10,13,14 IHE induced increased cross-sectional area (CSA) and increased coronary blood flow (CBF) in healthy subjects, but unchanged and decreased CSA and CBF were found in mild and severe CAD, respectively. 11 CEF-MRI measures were shown to be predominantly NO dependent, inversely related to underlying CAD, and reproducible. 14,15 Despite this significant advance, there are still limitations that hinder the widespread use of non-invasive CEF assessment by MRI. First, atherosclerosis is a diffuse and heterogeneous process, whereas 2D CEF-MR measures derive from single coronary segments and may miss the extent and exact locations of dysfunction. 11 Rather than two or three segments per study, coverage along proximal to mid-coronaries would allow for better CEF measurements along multiple segments. Second, current scans require 20-25 s breath-holds during exercise, which are difficult for some patients and limit spatial resolution (although sub-millimetric resolution is achieved in plane, the slice thickness remains about 8 mm). Third, the planning of 2D MRI for cross-sectional views of the coronaries requires expertise. Recent advances in self-gated whole-heart coronary magnetic resonance angiography (MRA) [16][17][18] and three-dimensional (3D) cine MRA [19][20][21][22][23][24] now allow for free-breathing acquisitions with isotropic resolution that do not require extensive planning expertise, although for high resolutions (<1.5 mm) scan time is either unpredictable or typically longer than the IHE duration ($7 min) of current CEF-MR studies.
To address the above limitations, we propose here 3D-Stars: a novel respiratory self-gated and resolved, volume targeted 3D coronary cine method to non-invasively assess CEF for the first time in 3D. 3D-Stars enables imaging of large proximal portions of the coronary artery, during free-breathing, with isotropic resolution, and in about 6 min-a time tolerated for handgrip stress.
3D-Stars uses a stack-of-stars 25,26 sequence modified for untriggered golden-angle-rotated 27 continuous imaging 19,20 and combines it with a cardiac and respiratory motion-resolved parallel imaging and compressed sensing reconstruction. 28 It includes a novel image-based respiratory self-gating approach, which is presented here and tested against a previously described method. 29,30 Additionally, we implemented a volumetric vessel analysis software tool that allows for semi-automated extraction of 3D CEF and image quality measures.
In this work, we test the ability of 3D-Stars to (A) image proximal and mid coronary arteries both at rest and during IHE stress by evaluating vessel length and sharpness and (B) derive non-invasive 3D CEF measures along the vessel in comparison with a reference 2D radial approach in healthy subjects.

| METHODS
The proposed 3D-Stars MRI approach represents an ensemble of MRI data acquisition, image reconstruction, and 3D vessel segmentation methods that allows for volumetric assessment of coronary anatomy and, potentially, vasoreactivity following provocation. First, a 3D stackof-stars acquisition volume is planned on localizers to image the proximal and mid-portions of the coronary artery. Data are acquired continuously and untriggered during free-breathing and later sorted into cardiac and respiratory phases using pulse oximeter gating and a novel image-based respiratory self-gating method, respectively. The sorted data are then reconstructed using five-dimensional (5D) golden-angle radial sparse parallel (GRASP) MRI 28,31 to reduce undersampling artifacts in the cardiac and respiratory motion-resolved images. From the resulting 5D dataset, a systolic or diastolic image can then be selected at end-expiration for volumetric analysis and display of the coronary artery, which includes the estimation of image quality parameters such as detected vessel length and sharpness as well as CSA along the centerline.

| Acquisition
3D-Stars acquisitions were performed on a 3 T clinical scanner (MAGNETOM Prisma fit , Siemens Healthcare, Erlangen, Germany) with a 30-channel cardiac coil array for signal reception. A targeted volume was prescribed on a whole-heart survey in a single oblique coronal orientation ( Figure 1A-C). The volume was placed to include the ostium ( Figure 1A) and all proximal to mid segments of the right coronary artery (RCA, Figure 1B and ) but excluded the left ventricle to ensure inflow of fully relaxed blood. 32 The prototype 3D-Stars sequence used a continuously acquired segmented hybrid radial (in-plane)/Cartesian (through-plane) gradient-echo acquisition. Each data segment comprised chemical shift selective fat saturation 33 followed by 20 radial views acquired along the slice direction (ie 20 slices) with center-out ordering and the same azimuthal angle to achieve fat-saturated contrast at the center of k-space. The duration of a segment was approximately 106 ms. Consecutive data segments were rotated by the golden angle (111.25 ) 27 on the azimuthal plane. Other protocol parameters included 1.44 mm isotropic resolution, 300 Â 300 Â 28.8 mm 3 field-of-view with single oblique coronal orientation, 3380 radial views/slice, 20 slices, T R /T E = 4.68/2.18 ms, 12 RF excitation angle and 361 s acquisition time independent of the heart rate. Figure 1D shows a schematic diagram of the acquisition.

| Respiratory self-gating
For the proposed image-based respiratory self-gating method, low-spatial/high-temporal-resolution 2D images were reconstructed retrospectively from the 3D-Stars dataset using compressed sensing from three adjacent central slices (voxel size = 7.5 Â 7.5 Â 4.2 mm) for each data segment with a 10-segment sliding window. This corresponded to a 106 ms frame update time and an approximately 1 s temporal footprint.
First, the 20 k-space lines of each segment were fast Fourier transformed in the slice-encoding direction (k z ) to select three lines from the three center slices within the volume. These three lines were then summed to obtain a four-dimensional (4D; k x -k y -t-c) time-resolved, multi-coil 2D k-space dataset, where k x , k y , t, and c represent in-plane k-space, time and coil dimensions, respectively. A sliding window of 10 consecutive data segments ($1 s) was applied along t. Supporting Figure S1 shows how self-gated respiratory image frames are generated. The central three slices of the 3D-Stars volume were used in all subjects to obtain a consistent view of the lung-diaphragm interface and sufficient signal for the reconstruction. This dataset was then fed to a k-t SPARSE SENSE (sensitivity encoding) 34 algorithm. Prior to reconstruction, data were compressed into five virtual channels 35 and coil sensitivity maps were generated with ESPIRiT (eigen-vector-based iterative self-consistent parallel imaging reconstruction) 36,37 from a central 8.5% portion of k x -k y using the first 1000 time points combined. The k-t SPARSE SENSE reconstruction included a multi-coil data consistency term and two total variation regularization terms (which enforce sparsity of finite differences) 38 in the image domain F I G U R E 1 Planning and acquisition of 3D-Stars. A-C, A targeted stack-of-stars volume was prescribed single oblique on a whole-heart coronary MR angiogram. Dashed lines indicate acquisition volume boundaries as they intercept transversal slices of the whole-heart survey. The volume was placed to include the ostium (A) and all proximal to mid segments of the RCA (B, C) but excluded the left ventricle to ensure inflow of fully relaxed blood. D, The stack-of-stars k-space was acquired continuously and untriggered for about 6 min. For each data segment, 20 radial views were acquired along the slice direction (ie 20 slices) with center-out line ordering and the same azimuthal angle, after chemical shift selective fat saturation. Consecutive data segments were rotated by the golden angle in the azimuthal plane and through time, respectively, and was solved in 40 iterations using an alternating direction method of multipliers (ADMM) optimizer. The total number of iterations was set empirically as the minimum number after which no changes were observed in the final reconstruction. Scaling factors for the regularization terms were found empirically and were both set to 0.01. These algorithms were implemented using the BART toolbox. 39 One-dimensional (1D) foot-head intensity profiles from a region of interest (ROI) drawn around the diaphragm were extracted from each image (Figure 2A, B). The displacement of the diaphragm during the acquisition was estimated with normalized cross-correlation, using a manually selected profile at end-expiration, after 16-fold Fourier interpolation for sub-millimeter resolution ($0.5 mm) ( Figure 2C). Thus, the estimated displacement reflects the relative respiratory position during the acquisition of each data segment. After discarding data segments at outlier respiratory positions (>2 standard deviations from the mean position), the remaining data segments were binned in eight or four respiratory phases with the same data size and without view-sharing ( Figure 2D).
For comparison, a previously described k-space-based respiratory self-gating method was also implemented. 40,41 This method extracts the k-space center amplitude of each data segment during the scan and recovers a respiratory signal from its fluctuation based on processing in both time and frequency domains. Further details are given in Supporting Figure S2. Note that this signal was only used for comparison in terms of correlation with signals from the proposed image-based respiratory self-gating method, and it was not used for data sorting or image reconstruction.
The processing for the image-based and k-space-based respiratory self-gating methods was implemented and retrospectively performed in MATLAB (MathWorks, Natick, Massachussets).

| Reconstruction of 3D-Stars
At the scanner, the entire 3D-Stars dataset was reconstructed and displayed as a single motion-combined (ie static) 3D image using vendor gridding reconstruction software modified to accept golden-angle rotated stack-of-stars data. These images were blurred due to motion but offered immediate feedback after the acquisition to check the planning of the 3D volume prior to IHE.
Offline, the full dataset was sorted along the time dimensions into 25 cardiac (t c ) and 8 respiratory (t r ) phases to obtain a 6D k-space dataset (k x -k y -k z -t c -t r -c) that was reconstructed with a 5D-GRASP algorithm 28,31 to yield cardiac and respiratory motion-resolved 3D cine images. For subjects with short quiescent periods during the heartbeat, 45 cardiac and 4 respiratory phases were used to better resolve cardiac frames with minimal coronary motion. The same data compression and coil sensitivity map estimation as described above were used. The 5D-GRASP algorithm included a multi-coil data consistency term and two total variation regularization terms 38 along the cardiac and respiratory time dimensions, respectively, and was solved in 20 iterations using a conjugate gradient optimizer. Spatial regularization was not performed. Scaling factors for the temporal regularization terms were found empirically and were set to 0.05 and 0.025 for cardiac and respiratory terms, respectively. This algorithm was implemented in MATLAB and based on source code available in Reference 42.

| Volumetric vessel analysis software
A semi-automated software for volumetric vessel analysis was implemented in the graphical programming interface (GPI) 43 to track the centerline and segment the lumen of the coronary in 3D. The tracked centerline provided a detected vessel length measure, allowed for estimation of CSA along the vessel, and also enabled multi-planar reformatting.
F I G U R E 2 A, B, The proposed image-based respiratory self-gating method used low-spatial/high-temporal-resolution 2D images reconstructed with a k-t SPARSE SENSE algorithm (see text) from three adjacent central slices (voxel size = 7.5 Â 7.5 Â 4.2 mm 3 ) for each data segment (every $106 ms) with a 10-segment sliding window ($1 s). C, 1D foot-head intensity profiles were obtained from an ROI including the interface between the lung and the diaphragm in each image (white rectangle) and were compared with normalized cross-correlation to estimate diaphragmatic displacement. D, Subsequently, this displacement signal was displayed over time to assign each data segment to a respiratory phase. Eight or four phases are created with the same data size For each 3D-Stars dataset of 25 Â 8 (or 45 Â 4) 3D images, one image showing the best depiction of the coronary is selected by visual inspection. Tracking of the vessel centerline was performed with a two-step method inspired by Soleimanifard et al 44 that accepts two user input start and end points along the coronary vessel ( Figure 3A). First, an initial centerline was estimated with a fast-marching method 45 from the user input points and a vessel-enhanced image as speed term. Here, a Frangi filtered image 46 was used as the speed term to produce an initial coarse estimate of the centerline. Then a 3D thresholding segmentation at 65% of the intensity range between a local maximum and a global minimum signal (both defined below) was performed along the current centerline to produce a coarse segmentation mask of the lumen, which was smoothed out using morphological operations. Second, a refined centerline estimate was obtained after another execution of the fast-marching method using same input points and a distance map of the current lumen mask as the speed term ( Figure 3B). The final centerline is obtained after smoothing with spline interpolation ( Figure 3C). Detected vessel length was defined as that obtained from such spline interpolation between the user input start and end point. If the software was not able to track the vessel between the two points, the end point was moved more proximal.
To estimate CSA along the vessel centerline, a 3D thresholding segmentation of the lumen was performed in adjacent 2 mm long tubes. The tubes were delimited every 2 mm by planes calculated to be orthogonal to the centerline. For the segmentation, thresholds of 50% and 65% of the intensity range between a local maximum and a global minimum intensity levels were tested. The local maximum was obtained as maximum intensity level in each tube with a 3 mm radius and 2 mm length. The global minimum was selected as the minimum image intensity in the 3D volume prescribed by the user start and end points. CSA values were estimated as (tube volume)/(tube length) from each tube along the centerline.
The 50%-threshold segmentation corresponds to a volumetric full width at half of the maximum (FWHM) method and was chosen in accordance with the 2D FWHM approach adopted in previous literature on CEF assessment by 2D MRI. 11,14,15,40,[47][48][49][50] The 65% segmentation, instead, corresponds to a volumetric full width at 65% of the maximum (FW65%M) method and was chosen to provide a more robust segmentation of the midportions of the vessel where the contrast and delineation of the lumen were found to be lower than in proximal segments in the 3D-Stars images.
A rendering of the entire FW65%M segmented volume in an example subject is also shown in Figure 3C. The volumetric CSA approach was validated in a numerical phantom, see Supporting Figure S3 for details.
Multiplanar reformatting was performed along the coronary ( Figure 3D), for soap-bubble visualization 51 and to evaluate vessel sharpness every 2 mm along the centerline, following Botnar et al. 52 Vessel sharpness values were averaged along the first 2 cm and in the remaining portion of the detected centerline length. Additionally, multiplanar reformatting perpendicular to the coronary centerline was performed to visualize crosssections ( Figure 3E).

| MR experiments
All human studies were approved by the Johns Hopkins School of Medicine Institutional Review Board and informed written consent was obtained from all study subjects. In vivo experiments were performed in eight healthy volunteers (four women, 29 ± 3 years, with no history of heart disease, hypertension or diabetes, and non-smokers).
3D-Stars CEF assessment of the proximal and mid portion of the right coronary arteries was performed in the morning after overnight fasting using the same IHE protocol as described previously, 11 and with a reference 2D radial cine technique for comparison. 49 F I G U R E 3 Main steps of the volumetric vessel analysis software for coronary lumen tracking, segmentation, and longitudinal and crosssectional reformatting. A, After the user has selected a 3D image frame in the 5D dataset with the best depiction of the coronary artery (usually in diastole and end-expiration), start and end points of the visible RCA are selected. B (left inset), From these two points, a modified fast-marching method tracks the vessel centerline. B (center and right insets), A volumetric full width at half (or 65%) of the maximum (FWHM or FW65%M) thresholding algorithm is used to segment the lumen. Within the segmented volume and along the centerline, adjacent 2 mm long cylinders are created to estimate the lumen CSA as cylinder-volume/2 mm length (displayed in alternating cyan and red color). C, A rendering of the entire segmented FW65%M volume and its centerline. D, E, The vessel can also be reformatted longitudinally for visualization and evaluation of sharpness (D), as well as orthogonally along the centerline (E) A whole-heart survey was first acquired to prescribe 3D-Stars with a single oblique orientation as mentioned above. On the same wholeheart scan, the three-point tool 53 was used to prescribe a targeted scan of the RCA, where 2D radial cine images were in turn prescribed perpendicular to straight segments of the proximal and mid-RCA. One free-breathing 3D-Stars acquisition and two single breath-hold 2D radial cine acquisitions (one at proximal and one at mid-level) were performed both at rest and during IHE. Thus for each subject, a total of two 3D-Stars acquisitions and four 2D radial cine acquisitions were performed. Each subject performed sustained IHE by continuously squeezing an MRI-compatible dynamometer (Stoelting, Wood Dale, Illinois) at 30% of his or her maximum grip strength under direct supervision. 12 Heart rate and blood pressure were measured at rest and during IHE stress using the scanner MRI-compatible ECG/pulse oximeter module and a noninvasive calf blood pressure monitor (Invivo, Precess, Orlando, Florida). The rate-pressure product (RPP) was calculated as systolic blood pressure Â heart rate.
The reference ECG-gated gradient echo 2D radial sequence was acquired as described before 49 with the following protocol parameters: spatial resolution = 0.9 Â 0.9 Â 8 mm 3 (to match gold standard CEF-MRI protocol 11 ), field of view = 260 mm, T R /T E = 4.9/2.7 ms, acquisition window (shot T R ) = 57.6 ms, 19 s breath-hold duration, spectral spatial water excitation, 22 RF excitation angle, 40 cardiac phases. Images were reconstructed at the scanner using the vendor gridding reconstruction.
Acquisition and reconstruction of 3D-Stars were performed as described in Section 2.1. During the free-breathing 3D-Stars acquisitions a pulse oximeter probe was applied to the toes (to free the hands for IHE) and timestamps were recorded for retrospective cardiac gating. Signal from a respiratory bellow was also recorded during the 3D-Stars acquisitions for retrospective comparison with image-based and k-space-based self-gating signals but was not used for data sorting and reconstruction.
All offline reconstructions, processing and analyses were performed using a Dell Precision T7810 workstation equipped with a 32-core Intel ® Xeon ® E5-2630 v3 CPU @2.40 GHz and 256 GB RAM.

| Respiratory self-gating
To assess the quality of the proposed image-based respiratory self-gating, the detected diaphragmatic displacement was compared with a previously described k-space-based self-gating signal and the signal obtained from the respiratory bellow. Before comparison, the respiratory bellow signal, as stored at the scanner with 200 Hz sampling frequency, was smoothed with a low-pass Butterworth filter of order = 6 and cutoff frequency = 4 Hz and down-sampled to match the self-gating signal temporal resolution (ie, acquisition time of one data segment ≈ 106 ms).
Respiratory gating signals were assessed using Pearson's correlation coefficient (CC) obtained by linear regression analysis between image-based and k-space-based self-gating, between image-based self-gating and respiratory bellow, and between k-space-based self-gating and respiratory bellow, with P < 0.05 considered statistically significant.

| 2D radial images
Image quality was assessed with a measure of vessel sharpness quantified at the interface of lumen and surroundings along 30 radial spokes centered on the lumen, following Botnar et al. 52 Vessel sharpness measures were performed after eightfold zero-filling in GPI and averages were determined across subjects for proximal and mid locations.
Coronary artery CSA measurements were performed using the semi-automated Cine tool (General Electric, Milwaukee, Wisconsin) based on the FWHM algorithm as previously described. 11 The percent change of CSA between rest and IHE was quantified as a measure of CEF. 11 Both vessel sharpness and CSA were determined in three diastolic frames, manually selected for minimal coronary motion, and averaged.

| 3D-Stars images
In the volumetric vessel analysis software, the selected image frame was interpolated by a factor of 8 with zero-filling in all three spatial dimensions prior to analysis. 48 Images from rest and IHE acquisition were analyzed separately, but the user start point was selected to match the same location, relative to the ostium of the RCA. Image quality was assessed with measures of detected vessel length, and vessel sharpness in proximal (first 2 cm of length) and mid-portions (>2 cm), as defined above.
CSA measures were obtained every 2 mm along the centerline with the volumetric FWHM and FW65%M segmentation methods. All CSA measures were smoothed using a moving average window of three samples (6 mm) and paired between rest and stress acquisitions at each location along the centerline. In the case of failure of segmentation, which can be due to insufficient contrast between the lumen and the surroundings or the presence of vessel side branches, leading to erroneous CSA estimate, both measures at same location from rest and stress acquisitions were excluded. Percentage change of CSA (%CSA-change) between rest and IHE was quantified as a measure of CEF. 11 Length of segmented lumen along the centerline was evaluated for both FWHM and FW65%M methods as a measure of segmentation performance.
To compare 3D-Stars with the reference 2D radial CEF measures, in this study acting as the gold standard, a location matched analysis was performed in which 3D-Stars FW65%M CSA and %CSA-change measures were selected at specific locations to match the proximal and mid-slice positions of 2D radial and compared with 2D radial FWHM measures. In order to compare 3D-Stars and 2D radial measures with the same FWHM segmentation threshold, a per subject analysis was performed, in which FWHM CSA measures were averaged in each subject for rest and stress acquisitions.

| Statistical analysis
Detected length (of 3D-Stars) and vessel sharpness (of 2D radial and 3D-Stars) were compared between rest and IHE stress acquisitions for each imaging method with paired two-tailed Student's t-tests. Additionally, vessel sharpness was compared between imaging methods for proximal and mid-portions with paired two-tailed Student's t-tests.
F I G U R E 4 A-C, An example of respiratory gating signals is shown from a handgrip stress acquisition. Circles indicate the respiratory positions of the acquired data segments; different colors indicate respiratory phase assigned by the self-gating method. The proposed image-based respiratory self-gating method measures the actual displacement of the interface between the diaphragm and the lungs in millimeters (A). This is beneficial, in comparison to k-space self-gating algorithms (B) or bellow signals (C), because it allows for the detection of abrupt changes of endexpiratory position (A, arrow). This is not the case for B and C, where data segments acquired after 250 s are assigned to the end-expiratory phases (cyan and magenta). D, A quantitative comparison shows good agreement of the proposed image-based self-gating with both k-space-based self-gating and bellow signals In the location matched analysis, CSA and %CSA-change measures of 3D-Stars and 2D radial were compared between rest and stress acquisitions with paired two-tailed Student's t-tests and linear regression plots. In the per subject analysis, FWHM CSA measures were compared between 3D-Stars and 2D radial with paired two-tailed Student's t-tests separately for rest and stress acquisitions.
For statistical tests, alpha = 0.05 was used as the significant threshold and a significant difference was indicated as *P < 0.05, **P < 0.01, ***P < 0.001. Values were reported as average ± standard deviation.

| RESULTS
Data were acquired for all 16 3D-Stars scans (one rest and one separate stress acquisition for all eight subjects) and all 32 2D radial scans (a proximal and a mid-location each acquired at rest and again during stress for all eight subjects). RPP measured at rest was 6456 ± 1294 mmHgÁbeats/s and increased by 24% during IHE stress (8007 ± 1102 mmHgÁbeats/s, p < 0.001).
For 3D-Stars, respiratory gating signals from image-based self-gating, k-space-based self-gating and respiratory bellow were successfully extracted in all 16 scans. Figure 4 shows an example where image-based self-gating allowed for measuring diaphragmatic displacement and assignment of respiratory phases even during an abrupt change in respiratory pattern ( Figure 4A, arrow), unlike the previous k-space-based approach 17,30,41 and the respiratory bellow. Overall, the proposed image-based self-gating method showed good agreement with the previous k-space-based approach and the respiratory bellow, CC = 0.83 ± 0.14 and 0.68 ± 0.15, respectively ( Figure 4D).
Example 3D-Stars images reconstructed at the scanner and offline are shown in Figure 5. In comparison to gridding ( Figure  phases, respectively. The acceleration factor was calculated by dividing the 6535 lines needed for a fully sampled image (=20 (slices) Â 208 (matrix) Â π/2 (radial sampling)) by 338 and 376 lines acquired per phase (=67 600 total lines of a 3D-Stars scan divided by 200 and 180 phases), respectively. The average respiratory bin widths are shown in Supporting Figure S4.
F I G U R E 5 Example 3D-Stars images as reconstructed at the scanner and offline. A, At the scanner, gridding reconstruction of unsorted and motion-combined dataset shows good depiction of static structures but blurring at the level of the heart. B, C, Images reconstructed after sorting in cardiac and respiratory phases with simple gridding and sum-of-squares coil combination to illustrate the amount of undersampling. D-F, In contrast, 5D-GRASP reconstruction shows improved image quality across motion phases and reduced streaking artifacts The selected 3D-Stars image frame for analysis corresponded to the coronary quiescent period during a diastolic phase and at or near endexpiration in all cases. Three out of eight stress acquisitions had low vessel sharpness (<20%) and thus were excluded from the image quality and CEF assessment (N = 5 remaining subjects were included). The success rate of 2D radial acquisition was higher, and only 4 acquisitions out of 32 could not be analyzed due to poor image quality.
Reformats of 3D-Stars diastolic frames at end-expiration showed good depiction of the RCA in longitudinal and cross-sectional views at both rest and stress ( Figure 6). Reformat of the volumetric FW65%M segmentation on the reformatted cross-sectional view shows good delineation of the lumen (green contour, Figure 6). 2D radial images of the same subject showed good image quality and sharp depiction of lumen cross-sections in both rest and stress conditions and for proximal ( Figure 7A) and mid-acquisitions ( Figure 7B).
F I G U R E 6 3D-Stars reformats of diastolic frames at end-expiration from two example subjects. Good delineation of the coronary can be observed on the longitudinal views both at rest and in stress. Insets show cross-sectional reformats at the level of the solid lines, and overlaid green contours reformatted from the segmented FW65%M volume show good delineation of the lumen F I G U R E 7 Images from reference 2D radial sequence acquired at proximal (A) and mid-(B) levels along the RCA (see "Planning" insets) in Subject 1 of Figure 6. Good image quality and sharp depiction of lumen cross-sections can be observed at both anatomical levels and in both rest and stress conditions. These views compare well with images obtained with 3D-Stars ( Figure 6A). A, The proximal cross-section was obtained 16.1 mm distal from the ostium as shown in the "Planning" inset, corresponding to 3 mm proximal from Figure 6A-b. B, The mid cross-section was obtained 28.8 mm distal from the ostium, corresponding to 3 mm distal from Figure 6A-d CSA and calculated %CSA-change measures along the coronary arteries with IHE, a measure of CEF, were obtained with volumetric segmentation of 3D-Stars images and are shown in Figure 9. We found that a 65% threshold successfully segmented the lumen along longer portions of the tracked centerlines compared to a 50% threshold, 26 ± 9.4 mm versus 13 ± 2.7 mm of total distance, respectively, during both rest and IHE stress, which corresponded to 13 ± 4.7 versus 6.5 ± 1.3 total segments per vessel for 65% versus 50% thresholds, respectively. For the analysis along the coronary (Figure 9) we therefore used CSA measures obtained with 65% threshold, and averaged them among all subjects where at least three measures from different subjects were available per location along the coronary. CSA trended higher during IHE than baseline ( Figure 9A), and the mean %CSA-change ranged between 5% and 30% ( Figure 9B), which agrees with literature values in healthy subjects during IHE 11,14,15,47 (Figure 9A and 9B). These CSA and %CSA-change average values were obtained in the proximal 27 mm of the coronary ( Figure 9C). In the Supporting Information, Figure S5 shows the length along the vessel of tracked centerline and segmented lumen (which provided CSA measures) for rest and stress acquisitions averaged across subjects.
For validation, CSA values obtained with 3D-Stars using FW65%M segmentation were compared with those from the standard 2D radial approach (using FWHM segmentation) at matched locations and are shown in Figure 10. A significant increase in CSA during IHE was detected with both techniques (Figure 10A and 10B), leading to CSA changes of 34.9% ± 19.3% and 16.7% ± 10.3% for 3D-Stars and the standard 2D radial approach, respectively ( Figure 10C). It should be noted that CSA values determined using 3D-Stars were significantly lower than those with 2D radial, and we suspected this was due to the higher segmentation threshold for 3D-Stars (65%) than for 2D radial (50%). When both 3D-Stars and 2D radial were analyzed with the same 50% threshold (FWHM), then there was good agreement between absolute CSA measures for the F I G U R E 8 Image quality evaluation on the five subjects with both rest and stress data. A, Detected vessel length of 3D-Stars obtained from the centerline tracking algorithm, from input start to end point, was not significantly different between rest and stress. B, Vessel sharpness of 3D-Stars averaged over proximal portions (first 2 cm) was higher than that in mid-portions (P < 0.01). However, it was not significantly different between rest and IHE stress within proximal or mid-portions. C, Sharpness for 2D radial images was not significantly different between rest and stress but higher in mid-than proximal portions (P < 0.001). 3D-Stars sharpness values averaged from rest and stress were similar to 2D radial in proximal portions but lower in mid-portions (P < 0.001). Bars indicate standard deviation two techniques ( Figure 10D). Supporting Figure S6 shows the correlations in CSA values and %CSA-change between 3D-Stars and 2D radial but the small sample size precluded rigorous statistical significance.

| DISCUSSION
3D-Stars is the first method to image the proximal and mid-portions of the RCA both at rest and during IHE stress, providing non-invasive assessment of CEF continously along the coronary artery with isotropic resolution. The isotropic and volumetric nature of the 3D-Stars datasets enabled quantification of CSA along several centimeters of the RCA both at rest and during stress. 3D-Stars is a free-breathing targeted volumetric coronary cine technique that combines a continuous stack-of-stars acquisition with a motion-resolved compressed sensing reconstruction exploiting cardiac gating and respiratory self-gating to obtain diastolic images at end-expiration of the RCA for quantitative analysis. A custom-built volumetric analysis software allows for semi-automated vessel centerline tracking and segmentation of the lumen to provide CSA along the centerline and F I G U R E 9 A, Paired CSA measures at rest and in stress from five healthy subjects could be obtained along 26 ± 9.4 mm of the coronary with FW65%M segmentation. CSA trended higher during IHE stress than at rest. B, The mean percentage change ranged between 3% and 30%, in agreement with the literature. Bars indicate standard deviation. C, Some measures at given distances could not be obtained due to failure of segmentation due to insufficient contrast or presence of a side vessel measures for detected vessel length and vessel sharpness. Local CEF measured by %CSA-change during IHE with 3D-Stars was in accordance with the literature and showed a trend similar to that of a conventional 2D radial approach. Vessel image sharpness measured with 3D-Stars was not significantly different from that measured with 2D radial in proximal portions, albeit lower in more distal, mid-artery segments.
In theory, the left anterior descending artery (LAD) could be studied with these techniques, but the orientation, inflow and other factors would require additional study and optimization. As such, this work focused on the RCA, with the highest SNR, and additional studies are needed to implement and optimize for LAD imaging.

| Acquisition
In order to quantify coronary artery lumen dimensions, images without motion blurring are required and thus it is important that image acquisition occurs during the period of the cardiac cycle with minimal coronary motion. Conventional ECG cardiac triggering allows one to do so by setting a fixed acquisition window duration with a fixed delay from the R-wave. IHE induced a 16% increase in heart rate in normal subjects on average, 11 and could cause a fixed trigger time set at baseline to acquire stress images at other than the quiescent period. To avoid the risk of missing the quiescent period, we here adopted a cine imaging scheme that allows accurate retrospective selection of a systolic or diastolic phase with minimal coronary motion, as has been applied in the 2D approach in numerous studies. 11,14,15,40,47,49,54 In this study, 8 out of 10 cases used a different cardiac phase during stress compared to rest with the 2D radial approach in the 5 subjects for proximal and medial RCA. The average heart rate in all 8 subjects scan sessions increased from 58 ± 9 bpm at rest to 67 ± 8 bpm during stress.
Additionally, the acquisition scheme fulfills the following requirements for a volumetric CEF assessment: isotropic resolution in 3D and an acquisition duration limited to handgrip time ($7 min). To fulfill all of the above, we chose a continuous golden-angle stack-of-stars acquisition that allowed for multi-phase volume-targeted imaging of the proximal and mid-coronary arteries. The stack-of-stars pattern was previously used for cine coronary imaging, 25,26 although without isotropic resolution, during breath-holding and with the use of contrast agents. Similar to the recent development in 5D coronary imaging, 19,20,31,55 the continuous data acquisition combined with golden-angle rotation of the interleaves F I G U R E 1 0 A, B, Comparison of CSA measurements obtained with reference 2D radial approach (A, using FWHM segmentation) and 3D-Stars (B, using FW65%M segmentation) at matched locations along the centerline. Increased CSA at stress was obtained with both techniques. Measures could be matched at 5 out of 11 locations for both rest and stress in three out of five included subjects. Lower CSA values measured using 3D-Stars compared with 2D radial could be due to the different threshold used for segmentation. C, Average %-CSA change from five matched locations trended higher for 3D-Stars. D, When measures were obtained with the same 50% threshold in both 2D radial and 3D-Stars, there was a good agreement in CSA between the two techniques takes advantage of the extra-dimensional (XD)-GRASP framework for reducing streaking artifacts and reconstructing a cardiac and respiratory motion-resolved dataset. However, while these studies adopt a whole-heart imaging strategy exploiting a 3D radial Koosh-ball-like trajectory, their imaging time is greater than 10 min and incompatible with tolerable handgrip time and many stressors. Here we adopted a stack-of-stars approach because of its shorter acquisition time due to significantly smaller FOV in the slice direction as compared with the Koosh ball acquisition. An additional 'pro' of this choice is the inflow contrast that can be taken advantage of. A drawback is that a scout scan and planning of the targeted volume are required and also that a reduced FOV may preclude visualization of distal portions of the vessel and tortuous side branches. 56,57 The need for fat saturation at the beginning of each data segment acquisition required a center-out line reordering scheme. This, in the adopted continuous acquisition, caused abrupt gradient changes between the acquisition of the last line of a given data segment and the first one of the next segment. Such gradient changes might be the cause for artifacts that we found on the ECG signal, preventing its use for cardiac synchronization. Because of these challenges encountered with ECG signals, we adopted a pulse oximeter for retrospective cardiac gating. As an alternative to fat saturation, water-excitation pulses could allow for different reordering schemes, since they no longer require acquisition of the center of k-space at the beginning of each data segment. For example, a line reordering scheme that alternates linear ascending and descending orders would limit gradient jumps in the slice direction to one step and likely reduce eddy currents and ECG signal artifacts. However, use of such water-excitation pulses or other fat insensitive techniques 58 was not pursued in this work.
The proposed method achieved a total scan time of about 6 min for a free-breathing 3D coronary cine with 1.4 mm isotropic resolution. This was chosen to exploit as much time as possible during IHE while allowing acquisition of 2D radial images for reference. Further development that shortens total scan time may reduce scan failure due to participant movement. The 1.4 mm isotropic resolution is worse than what is achieved inplane with the reference 2D radial approach and in the literature (0.89 mm) 11 . However, compared with the 2D radial approach (voxel size: 0.89 Â 0.89 Â 8 mm 3 ) 3D-Stars has a much smaller voxel volume (by $57%, [1.4 mm] 3 isotropic) and this may offer benefits, especially in the slice direction, which may in theory be helpful for the volumetric segmentation approach for estimating CSA.

| Respiratory self-gating
The proposed new image-based self-gating method for respiratory motion detection extracts 2D images from the 3D stack-of-stars acquisition about every 100 ms with an approximately 1 s sliding window. Considering a respiratory cycle duration of 4-5 s, the acquisition window length was found empirically as a trade-off between the sufficient temporal resolution to detect respiratory motion changes and the image quality that could be achieved with the described compressed sensing reconstruction, whose parameters were also optimized empirically to maximize contrast of the lung-diaphragm inteface. Therefore, the heart appeared blurred on these images, but was not used for motion detection; instead, the sharp lung-diaphragm interface was used to estimate respiratory motion states. The method accurately detected respiratory motion, as shown by the good agreement between the detected diaphragmatic displacement and the respiratory bellow signal ( Figure 4B) that was used for validation. The main advantage of this method is the ability to track the interface of the diaphragm with sufficiently high temporal resolution and measure its relative position in millimeters. This allowed for binning data from similar diaphragm positions, potentially leading to sharper images. We found two cases where the image-based approach detected a shift of the end expiratory position during the scan, whereas this was not detected with the bellow or the center of k-space method (one case shown in Figure 4A). Similar and strong correlation was found between each of the self-gating methods and the respiratory bellow, and a stronger agreement was found between the self-gating methods themselves ( Figure 4B). This suggests that both selfgating methods perform similarly well under regular breathing. An advantage of the implemented k-space-based self-gating method is the smaller footprint (one data segment) compared with the proposed image-based method, which uses a 10-segment sliding window. While 3D-Stars could be applied as a whole heart sequence, this extension would prevent the use of image-based self-gating in a protocol with many more slices, because its temporal footprint increases linearly with the segment length, and the k-space based self-gating method could be used instead. Another advantage of k-space-based self-gating methods is that they rely on band pass filters and analysis of the power density spectrum, which can usually be automatized. 40,41 In the proposed image-based method, the manual interaction required for drawing the diaphragmatic ROI for the image-based approach requires only a few seconds; however, this currently remains a limitation. 17,30,41 The method relies on the presence of the diaphram in the three central slices, and in our experience the dome of the diaphragm was clearly visible in all eight subjects. However, the efficacy of our imagebased method in more subjects and patients remains to be investigated. This work used the detected diaphragm position of each data segment in millimeters only for binning. However, this information could be used for both intra-and inter-bin motion correction to better align data segments, theoreticaly increasing sparsity between respiratory bins and improving the motion-resolved compressed sensing reconstruction as a result. 59

| Reconstruction
The motion-resolved reconstruction provided with the 5D-GRASP framework yields cardiac cine images at different respiratory positions.
However, the goal of this application is to select one image out of all cardiac and respiratory phases in diastole and end-expiration for quantitative CEF analysis. In four out of eight subjects, we used 25 cardiac and 8 respiratory phases, which allowed for adequate visualization of the diastolic period in two or three frames. This higher number of respiratory phases, as compared with previous work, 17 was adopted to have the option to choose between two frames for the end-expiratory phase in the case of drift or abrupt changes of the end-expiratory position. In the other subjects, a higher number of cardiac phases (fourty) was necessary in order to resolve the artery because of a shorter period of minimal coronary motion. For these reconstructions, we reduced the number of respiratory phases (to four) to keep similar data undersampling per frame and reconstruction computer memory load.
To enforce sparsity in 3D-Stars reconstruction, finite difference is used between neighboring temporal phases (temporal total variation) and reliably reduced undersampling artifacts ( Figure 5). Residual streaking was mostly confined to the outer FOV area and it did not affect centerline tracking to a noticeable extent. Other methods such as low-rank models 60 and patch-based reconstructions 61,62 could be investigated to exploit the high dimensionality of the 5D dataset and achieve a shorter acquisition time or increased FOV in the slice direction.
The use of finite differences in the spatial domain was avoided for the 5D-GRASP reconstruction of 3D-Stars images, as it is known to create a "patchy" look by reducing sharpness and flattening out signal due to a minimization process that penalizes finite differences in space.
Instead, such spatial total variation was adopted for the 2D respiratory self-gating images (in addition to temporal total variation), since these 2D images were used to visualize not small structures such as the coronary arteries but rather the marked lung-liver interface. We found that using finite differences in both space and time improved the 2D image quality while maintaining the sharpness of the diaphragm-lung interface.

| Analysis software and image quality
The analysis software tool provided for a user-friendly interface to select the best quality frame, to input start and end points along the coronary, and to perform all described post-processing steps. For each subject, the rest and stress datasets were displayed next to each other after eightfold zero-filling, and the start points for the segmentation were placed simultaneously on both datasets in order to visually match the beginning of the two RCAs. This procedure was important in order to compare the CSA measures between rest and stress at the same distance along the vessels.
Because of the particular FOV and orientation, it was difficult to find fiducial points at the beginning of the vessel; therefore, we opted for visually matching the two start points. Additionally, the start point has to be at a location where the lumen is cylindrical, therefore slightly distal from the ostium, for the tracking algorithm to work. To further minimize any errors in %CSA-change due to a potential mismatch between rest and stress, a 6 mm moving average of the CSA measures was applied along the vessel. Matching the locations of the start points on both rest and stress images could be improved in the future to minimize displacement between the two start points and make this procedure less operator dependent. A strategy to register the selected rest and stress datasets by placing one start point on an overlaid image of the two datasets could be explored in the future.
The vessel tracking algorithm was modified from its first implementation 44 by performing two iterations of the fast-marching method with different speed terms and allowed for a better centering of the centerline especially in distal portions. This, in turn, allowed for better determination of the orthogonal planes to create the 2 mm long tubes for analysis. With minimal user interaction of only two input points, the average vessel length measured by the algorithm was 60 ± 9.1 mm at rest and 49 ± 16.2 mm at IHE stress. These average values are lower in comparison with other studies 16,17,19,20,24,31,55,56,59,61,63 that reported visible vessel length above 80 mm but used manual selection of all points along the vessel. The vessels in our study were also visible beyond what the algorithm could track, and manual selection of all points 51 might have enabled measurement of more distal segments and a result closer to visual observations, but our goal was segmentation to estimate CSA. The ability of the fast-marching algorithm to track the vessel depends on the vesselness image (Frangi filter), which requires continuous contrast between the two selected points. When the algorithm failed to track the vessel, the end point was moved more proximal until tracking was successful. This is a limitation and future improvements of this process could include an automatic iterative search for the end point.
The vessel sharpness values measured in 3D-Stars were significantly lower in medial sections (distances > 20 mm) than in proximal portions. The fact that signal intensity of the lumen decreased with increasing distance along the vessel due to reduction of inflow effect could be a reason for such vessel sharpness and is a limitation of the current 3D-Stars method. On the other hand, 2D radial images leverage the inflow effect at every location, yielding high sharpness overall and significantly higher sharpness at the mid-level in comparison with 3D-Stars. However, 2D acquisitions do not offer volumetric coverage, require expert planning to obtain a slice perpendicular to the vessel, and the 8 mm slice thickness may be too large for tortuous segments. To account for this latter potential problem in 3D-Stars, CSA values were calculated from thresholding segmentation within 2 mm long cylindrical volumes along the centerline, and then planes perpendicular to the centerline were automatically created to separate the adjacent cylinders and guarantee correct segmentation even in tortuous sections. Albeit not significant, the sharpness of proximal RCA trended higher at rest than during stress, as can also be observed in Figure 6. In general, the long acquisition time (6 min) combined with the sustained handgrip exercise of stress scans might affect image quality.

| CEF assessment
3D-Stars provided measures of CSA for the first 27 mm at rest and during stress acquisitions in five subjects. Not all subjects contributed to each location, as indicated in Figure 9C, leading to non-uniformity of these measures in Figure 9A and 9B. While the %CSA-change average values are in the range reported previously for healthy volunteers, the standard deviations measured in this study were higher than previously reported. 11 Three subjects were excluded from the analysis due to low sharpness values of their 3D-Stars images. This could potentially be attributed to subject movement during the 6 min acquisition. For validation, 3D-Stars measures were compared with 2D radial acquisitions at matched locations.
Both techniques showed a significant increase of CSA with IHE stress, while CSA values for 3D-Stars ranged lower than those of 2D radial. One reason for these differences is the different thresholds applied to measure CSA between 3D-Stars and 2D radial. The 65% threshold was used in the analysis because it allowed for the detection of longer parts of the vessel, on average twice as long. FW65%M in 3D-Stars may be more robust in proximal to mid portions of the vessel because CNR is reduced as the inflow effect (and therefore signal, as described above) decreases distally. However, when the same threshold of 50% was used for both 3D-Stars and 2D radial, then CSA values were comparable between the two techniques ( Figure 10D).

| CONCLUSIONS
We developed and tested a novel coronary volumetric 3D cine MRI technique that exploits a golden-angle stack-of-stars trajectory for continuous untriggered acquisition, a novel image-based self-gating method for respiratory motion and a 5D-GRASP algorithm for motion-resolved reconstruction. "3D-Stars" allows one to obtain volumetric coronary cine data with 1.4 mm isotropic resolution, during free-breathing and in about 6 min. Additionally, a software tool was developed for the semi-automated analysis of coronary vessels that tracks the centerline and measures cross-sectional lumen area along the centerline with a volumetric thresholding segmentation. 3D-Stars was used to acquire non-invasive CEF measures along proximal and mid-portions of a coronary artery for the first time. In a small cohort of healthy subjects, 3D-Stars provided good image quality for the visualization of proximal and mid-portions of RCA at rest and during IHE stress, and yielded CEF measures along proximal portions of the RCA in accordance with the literature.