Volumetric Fetal Flow Imaging With Magnetic Resonance Imaging

Fetal development relies on a complex circulatory network. Accurate assessment of flow distribution is important for understanding pathologies and potential therapies. In this paper, we demonstrate a method for volumetric imaging of fetal flow with magnetic resonance imaging (MRI). Fetal MRI faces challenges: small vascular structures, unpredictable motion, and inadequate traditional cardiac gating methods. Here, orthogonal multislice stacks are acquired with accelerated multidimensional radial phase contrast (PC) MRI. Slices are reconstructed into flow sensitive time-series images with motion correction and image-based cardiac gating. They are then combined into a dynamic volume using slice-to-volume reconstruction (SVR) while resolving interslice spatiotemporal coregistration. Compared to prior methods, this approach achieves higher spatiotemporal resolution ( $1\times 1\times 1$ mm3, ~30 ms) with reduced scan time – important features for the quantification of flow through small fetal structures. Validation is demonstrated in adults by comparing SVR with 4D radial PCMRI (flow bias and limits of agreement: −1.1 ml/s and [−11.8 9.6] ml/s). Feasibility is demonstrated in late gestation fetuses by comparing SVR with 2D Cartesian PCMRI (flow bias and limits of agreement: −0.9 ml/min/kg and [−39.7 37.8] ml/min/kg). With SVR, we demonstrate complex flow pathways (such as parallel flow streams in the proximal inferior vena cava, preferential shunting of blood from the ductus venosus into the left atrium, and blood from the brain leaving the heart through the main pulmonary artery) for the first time in human fetal circulation. This method allows for comprehensive evaluation of the fetal circulation and enables future studies of fetal physiology.


I. INTRODUCTION
N ORMAL fetal development relies on an adequate supply of oxygen and nutrients from the placenta to critical fetal organs, such as the fetal myocardium and brain. Delivery is supported by a complex circulatory network that includes physiological shunts such as the ductus venosus, ductus arteriosus, and foramen ovale that are designed to both maximize oxygen transfer to the fetus and to protect oxygen delivery to the fetal brain. These shunts help deliver oxygenated blood from the placenta to the fetal brain and myocardium, with flows partially bypassing the fetal lungs and liver [1]. Fetal pathologies, such as fetal growth restriction and congenital heart disease, disrupt this distribution, causing injury to critical fetal organs and increasing risk of fetal mortality [2], [3]. Accurate measurement of this flow distribution is important because it may identify impairment of organs and thereby help to guide appropriate therapy and monitor efficacy. A complete assessment of the fetal circulatory system is possible through volumetric visualization, which can provide insight into the route of blood in the fetus.
Early studies demonstrated unique fetal flow mechanisms using injection of radiolabeled microspheres [4], [5]. Oxygenated nutrient-rich blood from the placenta partially bypasses the liver through the ductus venosus and flows into the proximal inferior vena cava (IVC p ). Upon arrival in the right atrium, this blood takes a preferential route through the foramen ovale into the left atrium, ultimately supplying oxygen and nutrients to the fetal brain and myocardium [5]. Two additional remarkable circulatory adaptations occur in this flow network. First, limited mixing of oxygen-rich blood from the ductus venosus with deoxygenated blood originating from the lower fetal body occurs in the IVC p resulting in two parallel streams. Second, deoxygenated blood from the lower and upper fetal body preferentially flows from the right atrium into the right ventricle with limited intracardiac mixing with blood from the ductus venosus. While these mechanisms support normal fetal development, pathologies that disrupt these mechanisms can impair growth and place the fetal brain at risk of hypoxic ischemic injury [2], [6], [7]. Assessment of flow distribution is hence crucial to monitor fetal growth and plan treatment.
Recently, phase contrast magnetic resonance imaging (PCMRI), whereby through-plane velocity is encoded in a single slice prescribed perpendicular to a target vessel, has been investigated for fetal blood flow quantification [8].
However, fetal PCMRI faces several challenges, such as small vasculature structures, high heart rates, uncontrollable motion, lack of traditional cardiac gating methods, and need for expert slice prescription during imaging. Recent advances in accelerated PCMRI have enabled quantitative assessment of the fetal circulation [2], [9]. These approaches have been used to quantify blood flow in great vessels, allowing estimation of fetal flow maps for healthy and pathological pregnancies [2], [9]. In cases of pathologies, redistribution of fetal flow has been tracked such that development of critical organs can be assessed and patient management can be better planned [2], [10]. However, these studies used single slice measurements with only through-plane velocity encoding. To achieve volumetric flow measurement, multidimensional velocity encoding (with each dimension being a velocity direction) is needed.
Conventional Cartesian 4D flow MRI, which involves timeresolved three-directional velocity encoding over an excited 3D volume, is well-established for creating flow maps through complex post-natal anatomy [11]. However, such methods are not well suited for fetal applications since they involve relatively long acquisition times during which gross or physiological fetal motion may corrupt the data. Gross fetal motions are unpredictable and vary in duration [12]. 3D reconstructions over short windows of consecutive acquisitions (real-time reconstructions) at sufficient spatial and temporal resolutions to resolve and correct such motions have not yet been demonstrated. Furthermore, characterizing the fetal heart rate from 3D acquisitions, for cardiac gating, is not possible because self-gating techniques are unreliable in the presence of fetal motion [13]. Alternatively, image based gating methods use real-time images at high temporal resolution to analyze dynamic features [14]. As explained later, the image-based fetal gating used in the current work relies on real-time images at high temporal resolution (∼25 ms) to resolve dynamic features. Finally, unrelated to motion correction, 2D imaging provides greater in-flow enhancement from blood flowing into the slice, which improves signal-to-noise and the contrast between blood and surrounding static tissue. This is important in fetal applications, which cannot receive endogenous contrast injection (a mainstay of postnatal 4D flow studies).
Multidimensional fetal PCMRI has been previously demonstrated using retrospective gating and motion correction [14] using golden-angle radial sampling [15]. With radial sampling, sufficiently high temporal resolution real-time images for motion correction and cardiac gating are achieved by updating the radial trajectory and the velocity encoding direction simultaneously, with each TR. While this technique provides high quality fetal images, it is constrained to single-slice use and provides only localized blood flow measurements. One solution to achieve volumetric fetal imaging in the presence of motion involves slice-to-volume reconstruction (SVR) [16]- [18]. The SVR approach, originally developed for fetal brain imaging and later adapted to fetal cardiovascular applications, acquires multiple overlapping stacks of slices which are combined using retrospective 3D motion compensation [12], [13]. Recently, 4D fetal whole heart CINE imaging has been demonstrated using SVR, along with the potential to provide velocity information in the final reconstruction [19], [20].
Here, we present a novel acquisition and reconstruction pipeline for volumetric fetal blood flow visualization and quantification at high spatiotemporal resolution using accelerated PCMRI. As explained in the following section, this pipeline relies on SVR to combine multislice radial PCMRI data acquired using an accelerated golden-angle acquisition [14]. This is the first demonstration of SVR using PCMRI, which provides many advantages over reported methods including an increase in spatial resolution by a factor of 6, an increase in temporal resolution by a factor of 2.4 and a reduction in scan time from 13 minutes to 7 minutes. This increase in spatiotemporal resolution with reduced scan time represents a major advance over prior methods for the evaluation of flow in small, complex anatomy. Furthermore, the demonstrated PCMRI method allows user-controlled flow sensitivity and provides intermediate 2D CINE reconstructions of flow, which can have clinical value independent of the final SVR reconstruction. Finally, this gradient-echo approach is less sensitive to off-resonance artifacts than methods based on balanced steady-state free precession (bSSFP). Validation of our novel approach is provided in adult volunteers, with feasibility demonstrated in four late gestation human fetuses.

A. Acquisition
Three orthogonal multislice stacks prescribed over a region of interest are acquired using multicoil accelerated radial multidimensional PCMRI, which is characterized by three orthogonal velocity measurements in a single slice. As described previously, it involves interleaved velocity encoding directions while continuously updating the trajectory angle by the golden angle [14]. In the stacks, slices have high in-plane resolution with relatively large slice thickness (see Methods for protocol details).

B. Reconstruction
Reconstruction involves four major steps, as detailed in the subsequent paragraphs and summarized in Fig. 1. First, each slice in each stack is reconstructed independently to yield fetal velocity sensitive cardiac-gated time-series images (CINEs) with retrospective motion correction and cardiac gating. Second, the slices are co-registered in space to correct for inter-slice motion. Third, the spatially co-registered slices are temporally synchronized such that the same cardiac phase in all slices occurs simultaneously in the reconstructed dynamic volume. Fourth, the spatiotemporally co-registered slices are reconstructed into a 4D fetal flow volume.

C. Step 1: Single Slice Reconstruction
In this step, each slice is reconstructed into a velocity sensitive CINE using a previously published framework [14]. The step for reconstructing one slice is explained below and the process is repeated for all slices in all stacks yielding velocity sensitive multislice CINEs. As depicted in Step 1 from Fig. 1, this involves real-time reconstructions for motion correction and cardiac gating followed by velocity sensitive Step 1: Reconstruction of slices with multidimensional phase contrast MRI encoding into 2D CINEs with motion correction, metric optimized gating, and compressed sensing (CS).
Step 2: Interslice motion correction. Stack-to-stack registration correcting for gross fetal motion with stacks of time averaged slices (#1) are registered to a reference stack (#2). Slices are first corrected for stack-to-stack registration (#3) and then iteratively corrected for finer inter-slice motions. A high resolution (HR) volume is estimated (#4) and slices are iteratively registered to it (#5). The registration parameters are used to update the HR volume (#6).
Step 3: Temporal registration of slices. Slices with largest overlap, based on the product of magnitude and speed, are found and registered in the time dimension while updating the temporally registered stack, L R , and temporally unregistered stack, L U . Step 4: 4D flow reconstruction by applying slice-to-volume reconstruction to the spatiotemporally co-registered flow sensitive slices. CINE reconstructions. Initial real-time reconstruction is performed for motion correction using a low temporal resolution (∼420 ms) window with compressed sensing (CS) [21]: where x moco is the low temporal resolution real-time reconstruction (moco denotes that the reconstruction is used for motion correction), F s is the coil-sensitivity weighted spatial Fourier operator, W moco selects temporal windows of spokes in data y, and p are sparsifying transforms with regularizing coefficients λ p . The sparsifying transforms used in this framework include spatial and temporal total variations (STV and TTV, respectively). Minimization is performed with a conjugate gradient descent algorithm [16].
To correct for in-plane fetal translational motion, real-time frames are co-registered using translational transformations over a user-defined region of interest spanning the fetal heart. Frames sharing low mutual information with other frames (less than the mean by 1.5× the interquartile range) are defined as through-plane motion and are automatically rejected from the pipeline. The remaining data represent quiescent periods, and the longest period is used for further reconstruction. Translational registration parameters are linearly interpolated to the repetition time (TR) of the acquisition and applied to correct motion in each spoke retained for reconstruction.
High temporal resolution real-time frames (∼ 25 ms) are then reconstructed using CS for cardiac gating: where x mog is the high temporal resolution real-time reconstruction, T moco is the motion correction operator, W mog creates temporal windows of spokes in acquired data y, and sparsifying transforms p are STV and TTV with regularizing coefficients λ p . The fetal heart rate is extracted using metric optimized gating (MOG) [22] followed by the data being sorted into cardiac phases. Finally, multidimensional velocity sensitive CINEs are reconstructed using CS: where c t is the multidimensional velocity sensitive CINE, G is a gating operator binning the data y into appropriate cardiac phase t, and sparsifying transforms p are STV and TTV with regularizing coefficients λ p . Velocity reconstructions are then corrected for background phase errors [23] and phase unwrapped for velocity aliasing.
The above steps are repeated for all slices in each stack, M i , to yield CINEs, c k,t i (which have an anatomical magnitude component m k,t i and 3 orthogonal velocity components, v k,t,d i with k the slice number, i the stack number, and d the velocity direction).

D. Step 2: Slice-to-Volume Registration
In this step, the reconstructed CINEs are time-averaged and co-registered spatially to each other. The time averaged stack with least inter-slice motion corruption is defined as the reference stack. This stack will have least data rejection, denoting least gross fetal motion, from the intra-slice motion correction step. Rigid registration is first performed to correct for gross motion between the stacks, yielding transforms which align each stack to the reference stack ( Fig. 1 Step 2).
A slice-to-volume registration is then performed iteratively to correct for the finer inter-slice motions. Rigid registration is performed to co-register each slice to a time-averaged volume which has been estimated from all the slices using SVR.
In the volumetric reconstruction, slices from different stacks are combined iteratively to generate an isotropic volume. The reconstruction attempts to invert the MRI sampling and motion processes which occurred during the acquisition. SVR is performed by minimizing the following cost function: where V I is the volume estimated using SVR, T k i is the rigid registration transform for the k th slice in the i th stack, R i is the stack-to-stack registration transform for the i th stack, N M is the number of acquired stacks, N i is the number of slices in the i th stack, B k i is the blurring operator based on the point spread function (PSF) of the acquisition process, and λ is the weighing coefficient for the 3-dimensional spatial total variation, STV 3D . The PSF of the radial MRI acquisition is jinc where dx, dy, and dz are spatial offsets from the reconstructed voxel, r is the in-plane resolution and σ z depends on the prescribed slice thickness. To make use of fast resampling functions, the in-plane PSF is approximated by a Gaussian kernel as in [16].
Slices are registered to V I to yield updated T k i rigid transforms. In this step, the registration solves for 3D translation and 3D rotation parameters which match each slice to its best location in V I . The latter are then used to update V I . The iterative process, registration followed by SVR, is repeated until the mutual information between the time-averaged volume and the acquired slices converges. In regions of incomplete slice coverage due to high fetal motion, interslice spatial interpolation is performed as part of the SVR. The final rigid transforms T k i at this stage provide spatial co-registration of all c k,t i slices, with corresponding velocity measurements v k,t,d i reoriented accordingly.

E. Step 3: Temporal Synchronization of CINEs
The spatially co-registered dynamic slices, c k,t i , must next be synchronized in time. To address this, temporal synchronization is performed across all slices based on their regions of overlap. Synchronization is performed using the product of image magnitude and speed, E k,t i given by The product of the time varying magnitude and speed gives a time varying angiogram. It is used for synchronization for the following reasons: (a) it leverages the correlated temporal fluctuations in both the magnitude and phase related to flow, and (b) phase noise in regions of low signal will be suppressed. Hence, this representation is robust to noise and accentuates the regions with blood flow for synchronization.
Strength of overlap between regions is first computed. An adjacency matrix for all slice pairs is populated using the sum of E k,t i in their region of overlap. A weighted graph is then computed representing the connections between slices [24]. The slice with the strongest connections is used as the seed point for the synchronization. As in Step 3 in Fig. 1, synchronization is performed iteratively by: (1) choosing a synchronized and unsynchronized slice with strongest overlap, (2) temporally shifting the unsynchronized slice until the entropy over the region of overlap is minimized, (3) looping until all slices are synchronized.

F. Step 4: 4D Velocity Reconstruction
With the slices, c k,t i , spatiotemporally co-registered, the slices are intensity corrected as in [17]. Scaling factors are computed between lowpass filtered slices in the reference stack and those from the other stacks. Intensity correction minimizes bias caused by anatomical magnitude disparities during combination of velocity sensitive complex data. The final 4D flow reconstruction (Step 4 in Fig. 1) is then performed by minimizing the cost function given by: The reconstructed volume V t I is a 4D volume containing dynamic volumetric magnitude and multidimensional velocity data. In this manuscript, 4D flow reconstruction generated using this approach will henceforth be referred to as 'radial SVR flow'.

A. Experimental Validation in Adults
To validate the framework, two adult volunteers (A1, A2; 25 years) were recruited with informed consent and scanned using a commercial 3T system (Prisma FIT , Siemens Healthineers, Erlangen, Germany) using body and spine matrices. Scanning protocol involved 3 orthogonal stacks of multislice 2D multidimensional radial PCMRI followed by a previously described 4D radial flow MRI acquisition [25]. Stacks of 2D PCMRI were performed in a sequential order and within each stack, slices were acquired in a sequential ascending order. This protocol was repeated 3 times in each adult. Scans were performed under free breathing conditions with the respiratory bellows and pulse oximeter gating physiological logs recorded. The prescribed volume focused on the main pulmonary artery, left pulmonary artery, right pulmonary artery, the aortic arch, and the descending aorta for radial SVR flow, and the whole chest for 4D radial flow. Relevant imaging parameters, chosen to keep scan times approximately equal between the two acquisition approaches, are summarized in Table I.

B. Fetal Scans
Four pregnant women with healthy singleton fetuses (gestation: 36 ± 1 weeks, subject identifiers: F1 -F4) were scanned with informed consent using the same commercial 3T MRI and multi-channel coil systems. For reference flow measurements, 2D Cartesian PCMRI acquisitions were performed on the major fetal vessels (DAO: descending aorta, AAO: ascending aorta, SVC: superior vena cava, MPA: main pulmonary artery, DA: ductus arteriosus) [22]. Acquisitions were repeated if an initial reconstruction on the scanner was deemed to be motion corrupted. Three orthogonal stacks of multislice multidimensional 2D radial PCMRI, centered on the fetal heart, were also acquired. Stacks of 2D PCMRI were performed in a sequential order and within each stack, slices were acquired in a sequential ascending order. Relevant imaging parameters are summarized in Table II. 500 radial spokes per flow encode were acquired, corresponding to an acceleration factor (R) of 12. Previous studies showed that radial acquisitions in human fetuses were accurately reconstructed with R ∼ 18 using CS [14], [26]. Hence, approximately 33% of the acquired data could be rejected during intra-slice motion correction while still allowing for CINE reconstruction.

C. Adult Flow Reconstructions
The acquired 4D radial flow data were first binned into 3 respiratory phases, based on the recorded respiratory log. The end expiratory phase was then binned into 15 cardiac phases, based on the pulse oximeter gating log, and reconstructed using CS (R ∼ 30) with spatial total variation (STV; λ = 0.008) and temporal total variation (TTV; λ = 0.01) regularization.
The multislice acquisitions were reconstructed as described above and in [14], but without intra-slice motion correction since acquisitions were short and initial real-time reconstructions showed negligible motion of the aorta and main pulmonary artery. Multidimensional velocity sensitive CINEs were then reconstructed using CS (R ∼ 13) with STV (λ = 0.001) and TTV (λ = 0.008) regularization with 15 iterations. Slices were then passed into the volumetric reconstruction pipeline, with inter-slice motion correction (with 1% tolerance in mutual information) and temporal alignment. Final velocity volumetric reconstructions were performed with STV (λ = 0.01) regularization for 10 iterations.

D. Fetal Flow Reconstructions
2D Cartesian PCMRI acquisitions were reconstructed into CINEs, with 15 cardiac phases, using MOG as previously described [8], [22]. The multislice radial acquisitions were reconstructed as described above and in [14]. Initial real-time reconstructions for motion correction were performed with CS using STV (λ = 0.01) and TTV (λ = 0.001) regularization for 15 iterations. Motion corrected data were reconstructed into high temporal resolution real-time images for cardiac gating with CS using STV (λ = 0.01) and TTV (λ = 0.001) regularization for 15 iterations. The fetal heart rate was extracted using MOG followed by binning the data into 15 cardiac phases (nominal temporal resolution approximately 30 ms, depending on fetal heart rate). Velocity sensitive CINEs were reconstructed with CS (R ≥ 12) using STV (λ = 0.01) and TTV (λ = 0.05) regularization for 15 iterations. Slices were then passed into the volumetric reconstruction pipeline for inter-slice motion correction (with 1% tolerance in mutual information) and temporal alignment. Final velocity volumetric reconstructions were performed using STV (λ = 0.01) regularization for 10 iterations. In regions of incomplete slice coverage due to high fetal motion, spatial interpolation was performed as part of the SVR.

E. Adult Analysis
All flow analysis and visualizations were performed using a prototype software (4D Flow v2.4, Siemens) [27]. The heart and great vessels were first segmented by thresholding the 4D flow data based on peak velocities, to create a phase-contrast angiogram. Refined semi-automated segmentation was then initiated through user-selected seed points for vessel centerline and lumen extraction. Mean flows were recorded at similar positions along the aorta, main pulmonary artery, left pulmonary artery (LPA), and right pulmonary artery (RPA) in the 4D flow and the radial SVR flow reconstructions. The flow values were compared using a paired t-test (with significance level set at p < 0.05), linear regression, and Bland-Altman analysis (with the difference computed by subtracting the radial SVR flow from the 4D flow measurements). The same analysis was repeated on only the smaller vessels (LPA and RPA) which are representative of SVR for the fetal vessels.

F. Fetal Analysis
2D Cartesian reconstructions were analyzed using the commercially available software QFlow ® from Medis Suite MR (Medis Medical Imaging, Leiden, The Netherlands).
From the reconstruction pipeline for the radial MRI data, the mean and range of the displacements and the percentage of data rejected owing to through-plane motion during intra-slice registration were recorded. The mean and range of the displacements and rotational Euler angles were also tracked for inter-slice registration. The mean and the standard deviation of the fetal heart rate obtained from MOG were recorded.
All volumetric flow analysis and visualizations of radial SVR flow data were performed using the prototype software as before. For flow measurements, contours were initially placed at 3 locations along each visible great vessel. The contours were refined by reorienting them perpendicular to the vessel axis and the flows through the cross sections were measured. The mean flows, indexed to fetal mass [8], from the volumetric and Cartesian reconstructions were compared by paired t-test, linear regression, and Bland-Altman analysis (with the difference computed from Cartesian PCMRI measurements subtracted from radial SVR flow). Internal consistency checks in the volumetric reconstruction were performed by assessing the conservation of mass in the reconstructed circulatory network. The flow in the DAO was compared with the incoming flows from the DA and the aortic isthmus (AI) through DAO = AI + DA. The flow in the SVC was compared with the flow leaving the aortic arch through SVC = AAO -AI. The flow entering the heart from the SVC, IVC (which is assumed to be equal to DAO flow), and the pulmonary veins (PV, given by the difference between the MPA and DA flows) were compared with the combined ventricular output (CVO = 1.03 × [MPA + AAO]) [1]. The consistencies were computed as a percentage of the CVO.

G. Fetal Flow Visualization
Radial SVR flow allowed for various ways to visualize the fetal circulation. Representative examples, depicting volumetric flow distribution, from the four fetal cases were reported here. To visualize fetal blood flow, particles were seeded throughout the segmentations in subject F1 and F2. Path lines were generated to visualize the movement of blood in the fetal heart and surrounding great vessels. To track fetal blood flow originating from specific locations, particle tracing was performed by placing emitter planes in target fetal vessels. In subject F3, emitter planes were placed in the ductus venosus and IVC to visualize the paths of oxygenated blood coming from placenta and deoxygenated blood coming from the lower fetal body. In subject F4, emitter planes were placed in the SVC and IVC p to look at the path of blood in the right side of the heart as it arrives from the fetal upper and lower body.

H. Reconstruction Specifications
All reconstructions in this work were performed on a standard desktop computer with the following specifications: 32 GB memory, processor Intel Core i7-6700 (3.40 GHz, 8 cores) and GPU Nvidia Geforce GTX 960 (2GB and 1024 CUDA cores). Regularization coefficients (λ) were optimized through empirical testing. For 4D radial flow, values for λ were chosen such that reconstructed flow values were internally consistent. Values of λ for CS were chosen such that the undersampled reconstructions were consistent with highly sampled reconstructions in prior tests. Since the acquisition of highly sampled dataset for SVR is impractical in-vivo, λ for SVR was solved using a grid search through simulated data by finding the λ that yielded a reconstructed 3D phantom with lowest errors from simulated slices.

A. Adult Flow Comparisons
4D radial flow reconstructions in adults required approximately 3 hours. Radial SVR flow reconstructions required approximately 5 hours.
Representative reconstructed CINE frames from one subject (A1) are presented in Fig. 2A, showing anatomy and flow in peak systole from each orthogonal stack. Subsequent SVR reconstruction based on multislice CINE data found that inter-slice displacements and rotations were small across all  Fig. 2B shows a rendering of the 3D cardiac anatomy from A1 reconstructed with SVR. Oblique resampled slices, from the reconstructed volume, depicting the surrounding great vessels near the heart are shown in Fig. 2C. Particle traces tracking blood flow in the aorta and MPA evolving over a systolic phase are shown in Fig. 2D. Blood can be seen going into the left pulmonary artery from the MPA and down the DAO from the AAO. Comparison of mean flows between the radial SVR flow and 4D radial flow reconstructions is shown in Fig. 3. Across all vessels, linear regression (Fig. 3A)

B. Fetal Flow Comparisons
Reconstruction of each slice from the multidimensional radial phase contrast MRI in the fetus required around 45 mins  slice are depicted in Fig. 4A. The differences between the extracted heartbeat durations obtained with MOG in presence and absence of motion correction are shown in Fig. 4B. Differences in the reconstructed CINEs at systole and diastole in the presence and absence of motion correction are presented in Fig. 4C. The septum is more conspicuous and cardiac features are less blurred with motion correction in both cardiac phases. The cardiac cycle duration across all fetal subjects was 430 ± 25 ms. Table III summarizes the intra-slice motion estimates, amount of data rejection and heart rate for the fetal data calculated by the reconstruction pipeline.
Representative reconstructed CINEs from one fetus (F3) are shown in Fig. 5A. The anatomy and velocity reconstructions at systole are presented for fetal axial, coronal, and sagittal planes. From SVR based on the CINE data, the mean interslice displacement was 5.6 ± 4.6 mm for all slices across all subjects. The corresponding mean rotational Euler angles were E x = −0.4 • ± 0.9 • , E y = −0.4 • ± 0.6 • , E z = −0.1 • ± 0.7 • . Largest slice-to-volume displacements were observed in F3 (7.7 ± 6.8 mm, maximum: 18 mm, angles [minimum maximum]: Fig. 5B shows a rendering of the 3D cardiac anatomy and surrounding great vessels in the F3 reconstructed with SVR. Oblique resampled slices, from the reconstructed volume, are shown in Fig. 5C. One visualization result of the SVR is shown in Fig. 5D using particle traces to track the ejection of blood from the ventricles into the aorta and MPA over a systolic phase. Blood can be seen going from the right ventricle into the MPA and from the left ventricle into the AAO. For comparison, the anatomical and through-plane flow encoded reconstructions from the 2D Cartesian PCMRI acquisition for the MPA are shown in Fig. 5E. The improvements with the inter-slice motion correction step in SVR are illustrated in Fig. 6. The improvements brought by each iteration of the interslice registration are shown in Fig. 6A for two representative planes through the reconstructed volume. For instance, the septum is becoming sharper and more conspicuous. The cardiac chambers are becoming more well defined. Vessels are becoming more conspicuous. The effects of motion correction on segmentation and flow analysis are shown in Fig. 6B. Combination of slices in SVR without inter-slice motion correction shows different anatomical regions being merged. Sections of the anatomy were blurred or missing from this reconstruction relative to the segmentation from the reconstruction with motion correction. Without inter-slice motion correction, the ventricles are smaller and have rougher edges with the main pulmonary artery unconnected to the right ventricle. The superior vena cava is not trackable in the absence of motion correction. Sections of the ascending aorta are not visible for segmentation. Table III summarizes the inter-slice motion estimates for the fetal data calculated by the reconstruction pipeline Fetal flow waveforms obtained using radial SVR flow showed good agreement with 2D Cartesian PCMRI measurements, as demonstrated by the representative curves from one fetus (F3) shown in Fig. 7. Features of these waveforms follow typical fetal arterial and venous patterns, and are consistent with prior publications: for instance the biphasic pulsatile behavior in the SVC can be observed [8], [28].  Fetal mean flows measured using radial SVR flow compared well with those obtained from target vessels using 2D Cartesian PCMRI. As shown in Fig. 8A, linear regression between flows across all vessels and fetuses had a slope = 0.94 (95% confidence interval = [0.72 1.16]), intercept = 8.9 (95% confidence interval = [−29.3 47.2]), and r 2 = 0.86. Bland-Altman analysis (Fig. 8B) had a bias of −0.9 ml/min/kg and limits of agreement [−39.7 37.8] ml/min/kg. Paired t-test showed that there was no significant difference between the measurements ( p = 0.86). The CVO across the subjects was 385 ± 41 ml/min/kg from the radial SVR. The corresponding measurement from the 2D Cartesian PCMRI was 375 ± 46 ml/min/kg. Table IV compares the mean flows and ranges of flows in fetal great vessels from literature [9], with Cartesian PCMRI and radial SVR flow. The measured flows in all the categories fall in the same range. Conservation of mass checks across the various cardiac structures showed differences of 7 ± 9%. Proximal to the aortic isthmus, the percentage difference between DAO and AI + DA flows, relative to CVO, was 8 ± 5%. The percentage difference between the SVC and the flow leaving the aorta for the upper body was 11 ± 8%. The difference between the CVO and the blood entering the heart was 2 ± 14%.

C. Fetal Flow Visualization
Particle traces tracking fetal blood flow in the heart and outflow tracts are depicted in Fig. 9A for F1 and Fig. 9B for F2. Ejection of the blood from the right ventricle into the MPA and the left ventricle into the AAO are depicted in Fig. 9A1 and Fig. 9B1. Blood flow from the SVC into the heart can also be observed in Fig. 9A1. Blood moving into the ventricle during diastole can be observed in Fig. 9A2 and Fig. 9B2. Supplementary File 1 shows a movie for the flow in F2 at different cardiac phases.
Particle traces depicting flow in the proximal IVC are shown in Fig. 10A with planes set in the ductus venosus (red) and distal IVC (blue) simultaneously emitting particles in F3. The figure demonstrates the path of blood over two heart beats. While both streams meet in the proximal IVC, the blood streams from the two vessels have limited mixing (Fig. 10A1). Blood from the ductus venosus preferentially shunts to the left side of the heart and most of the blood from the distal IVC goes to the right side of the heart (Fig. 10A2). Blood that originated from the ductus venosus or distal IVC exits the heart through the AAO or MPA respectively (Fig. 10A3). Fig. 10B depicts particle traces of flow from SVC and proximal IVC in F4. Simultaneous emission of particles from planes in the SVC (blue) and proximal IVC (red) are shown in Fig. 10B1. SVC blood flows into the right side of the heart while that from the IVC flows into both sides of the heart (Fig. 10B2). Traces in the AAO show that blood originated only from the IVC while those in the MPA show that the flowing blood originated from both the SVC and IVC. Supplementary File 2 shows a movie for the tracked flow paths in F3 and F4 at different cardiac phases.

V. DISCUSSION
In this study, we devised an approach using rapid multislice multidimensional radial PCMRI to image the volumetric human fetal blood flow. Reconstructions were performed using CS with cardiac gating through MOG, and retrospective intra-slice and inter-slice motion correction. Radial SVR flow was tested in adult subjects and compared with 4D radial flow MRI measurements. Feasibility of the approach was demonstrated in four human fetuses in the late third trimester and compared with 2D Cartesian PCMRI measurements in the great vessels. Along with fetal blood flow quantification, we were able to visualize complex fetal hemodynamics. We showed parallel streams of blood in the IVC and preferential blood flow routes in the fetal heart.
In the experiments in adults, radial SVR flow showed good agreement with 4D radial flow reconstructions. Since radial SVR flow was based on 2D slice acquisitions, it benefited from in-flow related enhancement with saturated spins being washed out of the imaged region and fresh spins moving in. This allowed the blood vessels to have higher intensities and become more conspicuous in the SVRs.
To achieve good contrast in conventional clinical 4D flow MRI, contrast agents (such as ferumoxytol and gadolinium) are injected [29]. However, contrast agents have disadvantages. In addition to their expense, injection of the contrast agent into the fetal circulation is challenging owing to potential injuries in case of fetal motion and difficulty in localizing the site of injection. Moreover, ferumoxytol can lead to hypertension and gadolinium has been linked to nephrogenic systemic fibrosis in patients with kidney failure [30], [31]. With such drawbacks, contrast agents cannot be administered to fetal subjects and the flow enhancements in the SVR approach provide an appropriate tool for imaging the in utero circulatory system.
In this work, we imaged the human fetal circulation with MRI at higher spatial and temporal resolutions than previously demonstrated, using a novel PCMRI acquisition and reconstruction approach. High spatial resolution is needed to resolve flow through small structures in the fetal anatomy, and this becomes increasingly important at earlier gestations and in complex congenital heart disease. Similarly, high temporal resolution is needed to resolve dynamic flow at fetal heart rates and consequently to trace the path taken by blood, as recommended for 4D flow measurements in postnatal applications [32]. We demonstrate the feasibility of using this approach to quantify flow in human fetuses at 1 mm isotropic spatial resolution and approximately 30 ms temporal resolution. To achieve such high resolutions within a clinically acceptable imaging time of approximately 6 minutes, data were acquired at high acceleration factors and reconstructed using CS. Previously reported methods for human fetal 4D flow have relied on data acquired using 2D k-t SENSE real-time bSSFP, achieving final reconstructions at 2 mm isotropic spatial resolution and 70 ms temporal resolution [20]. Increasing these resolutions may be possible using a higher k-t factor, with the corresponding risk of increased aliasing and reconstruction artifact. As with the method proposed here, higher resolutions will produce lower SNR, making real-time imaging and motion correction more challenging. The greater sensitivity of bSSFP to off-resonance, versus the gradientecho PCMRI method used here, may also pose challenges, especially at higher magnetic field strengths. Using multidimensional PCMRI for imaging the fetal circulation provides two advantages over previous literature. Firstly, the velocity sensitivity of the PCMRI approach may be controlled by the user. This is more challenging for methods that rely on the incidental flow sensitivity of other imaging sequences like bSSFP [33]. Secondly, our intermediate reconstructions before SVR provide flow sensitive 2D CINEs which have independent clinical value.
The range of flow values observed in this study using radial SVR flow agreed well to that obtained using Cartesian PCMRI and reference normal fetal blood flow measurements reported in literature [8], [9]. Differences between the radial SVR and Cartesian PCMRI flow curves are owing to two factors. First, there were time delays up to 30 minutes between the Cartesian PCMRI acquisitions and the radial multislice acquisitions. Physiological variations over these periods can lead to differences in measured flows. Second, the temporal resolution of the Cartesian PCMRI acquisitions was lower than the radial SVR acquisitions by a factor of two. Lower temporal resolution can lead to underestimation of peak flows. Furthermore, the values for internal consistency based on conservation of mass were comparable to internal consistency checks performed with 4D flow in animal fetuses [34]- [36]. These observations build confidence in the conservation of mass checks in this study. Also, it is important to note that the Cartesian scans were repeatedly acquired until an initial reconstruction on the scanner did not show signs of motion corruption. On average, Cartesian flow measurements were acquired twice, with localizer scans also rerun for prescription. Consequently, the total imaging time of the Cartesian scans (∼ 10 mins) was longer than the imaging time of the radial SVR flow (∼ 6.5 mins) scans.
In addition to the quantitative agreements, the volumetric visualizations corroborate previously demonstrated circulatory physiological mechanisms in animal fetuses [5], [34]. We provided visualization of parallel streams in the IVC in the human fetus, including preferential streaming of oxygenated blood from the ductus venosus through the foramen ovale and deoxygenated blood restricted to the right side of the heart. Moreover, radial SVR flow also showed that deoxygenated blood returning from the brain through the SVC is restricted to the right side of the heart and returned to the placenta via the MPA, DA, and DAO. Overall, these observations demonstrated the complex mechanisms in the fetal circulation that shunts oxygenated blood preferentially to the developing brain and heart. These findings support the idea that 3D visualization of the fetal circulation provides a more complete assessment of fetal circulation and its development than with 2D PCMRI. This technique provides opportunities to study a range of diseases that can place the fetal brain at risk of ischemichypoxic injury. For example, it may aid our understanding of various forms of congenital heart disease that impair brain development by altering the streaming of oxygenated blood flow to the fetal brain. Furthermore, by correlating MRI-derived observations of progressively fragmented fetal streaming, clinicians may be able to better plan treatments in-utero and postnatally.
Despite the success of the pipeline, there are limitations to this study. First, fetal radial SVR flow scans were part of a larger clinical research study. As a result, scan time allocated for this component was limited to approximately 6 minutes. This allowed acquisition of only three multislice stacks and without slice overlap, which is less than that used and recommended in similar studies of the fetal brain using SVR reconstruction [16], [19], [37]. This limited the slice coverage. Future work will acquire overlapping slices with larger anatomical coverage to remove the need for inter-slice interpolation. Second, SVR used in our approach depends on reconstructed CINEs which in turn depend on the amount of data rejected for intra-slice motion correction. Compared to other fetal volumetric reconstruction techniques, this limits our approach in cases where fetal motion is too severe and large portions of data must be rejected (for instance, 33% of data can be rejected to reconstruct CINEs [14]). Third, golden-angle radial sampling may suffer from data clustering once data are sorted into cardiac phases [38]. Clustering can result in large unsampled regions in k-space which lead to severe artifacts in the reconstructed images. In future implementations, different sampling patterns will be investigated to counteract this issue.

VI. CONCLUSION
In conclusion, we have developed an approach for dynamic volumetric fetal flow imaging using highly accelerated multidimensional radial PCMRI. The method makes use of real-time reconstructions for intra-slice motion correction and derivation of fetal heart rate to provide flow sensitive CINEs. They are then combined with inter-slice spatiotemporal registration and volumetric update to create 4D flow reconstructions in human fetuses. This work resulted in the visualizations of complex fetal cardiovascular physiology, such as parallel streaming in the proximal IVC and preferential routes taken by oxygenated blood in the fetal heart. It provides an important means of assessing fetal circulatory system and may enable new ways to manage fetal pathologies.