2.1.

Data Acquisition

The BOLD fMRI data were acquired from eight pregnant women with healthy fetuses and eight pregnant women with fetuses diagnosed with CHD between 25 and 40 weeks of gestation ($33.29\pm 4.08$). The study was performed by the Department of Radiology at the Children’s National Medical Center (CNMC) and approved by the CNMC Institutional Review Board. Written informed consents were obtained from all study participants.

Echo planar imaging (EPI) sequences were acquired on a 1.5T MR scanner (GE Healthcare, Inc., Waukesha, Wisconsin) using an eight-channel receive-only cardiac coil (USAI, Inc., Aurora, Ohio) with the following parameters: a matrix size of $128\times 128$, in-plane field of view (FOV) $420\times 420{\mathrm{mm}}^2$ or $440\times 440{\mathrm{mm}}^2$, slice thickness 5 to 8 mm with between-slice gap of 2 mm, repetition time (TR) of 2000 or 3000 ms, echo time of 1000 ms, a flip angle of 90 deg, and 144 or 288 volumes per acquisition. Each EPI sequence was obtained in an interleaved slice order to minimize cross-talk between adjacent slices. The fMRI acquisition protocol included both the placenta and fetal brain as shown in Fig. 1. The data were acquired in both tilted-coronal and axial planes to assess the effect of imaging plane on motion correction. The protocol parameters for both acquisition planes are summarized in Table 1.

Fig. 1

An example of fMRI data acquired from a healthy fetus in the coronal plane. It includes both the placenta (in red) and the fetal brain (in green).

Table 1

Protocol parameters for both axial and coronal acquisition planes.

The hyperoxia task paradigm consisted of three phases: a 2-min normoxia (baseline, 21% oxygen), followed by 100% oxygen ($15\mathrm{L}/\min$) administered via a maternal facial oxygen mask for 3 to 4 min, and an additional 2:12 to 3:36 min of normoxia (21% oxygen) to quantify return to baseline.

2.2.

Preprocessing Pipeline

As shown in Fig. 2, the fMRI data were processed according to the following pipeline which consists of four steps: (1) bias field correction, (2) motion correction, (3) VOR, and (4) data imputation. First, the inhomogeneous distribution of MR receiver sensitivity was corrected using the four-dimensional nonparametric bias estimator based on the B-spline field approximation^32,33 previously used for BOLD signal analysis of in vivo fetal fMRI of the brain.¹²

Fig. 2

The proposed preprocessing pipeline for slowly cycled block design fMRI of the fetus.

The fMRI dataset was temporally divided into three independent phases (baseline, hyperoxia, and return to baseline) according to the block design paradigm. The initial five volumes were excluded from each phase, since these are usually affected by phase transition. The phase-specific global motion correction (GMC) was carried out separately in each phase using six degrees of freedom rigid-body image registration using FLIRT as a part of the FMRIB Software Library (FSL).³⁴ The registration parameters were as follows: normalized cross-correlation ratio cost function, 256 histogram bins, trilinear interpolation, and search angle of $-20\mathrm{deg}$ to 20 deg in $X$, $Y$, and $Z$ axes. The temporal mean ($t$-mean) volume of each phase was used as the reference or template for image registration.

For local motion correction (LMC), two ROI masks of the placenta and fetal brain were manually defined from the $t$-mean volume of each phase using ITK-SNAP.³⁵ A total of six ROI masks of the brain and placenta (two for baseline, two for hyperoxia, and two for return to baseline) were created. The quality of the created masks was verified and corrected by an experienced clinician scientist. The ROI masks were dilated using a $3\times 3\times 3$ box kernel (in voxels) to allow a tolerance for the search window of local motion into a neighborhood of the ROI. In each phase-specific fMRI sequence, the FOV was restricted to the ROI masks of the fetal brain and placenta. LMC was then performed in each ROI. The LMC parameters were the same as the ones used in GMC. Both global and LMC steps were accelerated using sun grid engine parallel processing based on graphic processing units of the supercomputing system so that image registration could be performed simultaneously on multiple volumes.

The outlier volumes, which were significantly misaligned even after motion correction, were automatically detected and rejected by setting all voxels of those volumes as zero. To determine if a volume was rejected or not, the temporal outlier score (TOS) was computed for all voxels in the ROIs using the program 3dToutcount from AFNI where any signal falling outside 1.5 times the interquartile range (IQR) was regarded as an outlier. The motion outlier probability (MOP) was then computed by counting the number of outlier voxels, and the volume outliers were determined by thresholding the MOP values. Following VOR, all three phases were concatenated for each ROI again, and the BOLD fMRI signals were averaged as the median value over all voxels in each ROI. In this step, outlier voxels with abnormal intensities due to either physiological noise or other MR artifacts were automatically detected and excluded from the ROI-averaging process. Finally, all the missing data (set to be zero by VOR) in these ROI-averaged time series were replaced with the values estimated through statistical data imputation based on local polynomial regression whose kernel bandwidth ($s$) and the degree of polynomials ($p$) were set to be $s=40$ and $p=3$, respectively. The theoretical details underlying the above preprocessing pipeline are summarized in the next section.

2.3.

Design-Optimized Motion Correction

The traditional approach to motion correction in fMRI data of the brain is to align a volume to a predefined template using a geometrical transformation.^34,36,37 Without imposing any constraints, it can be treated as an image registration problem of finding the optimal geometric transformation that maximizes the similarity between volume and template. However, this typical approach is not directly applicable to fMRI data that include independently moving objects (i.e., the placenta and fetal brain) and are acquired with the epoch-based paradigm of sequential block design (i.e., short-term maternal hyperoxia) where the frequency of phase transition is slow.

To address such practical limitations of these unique motion correction challenges in stimulus-based fMRI of the moving fetus, we propose the design-optimized motion correction (DOMC) framework that is specific to the experimental paradigm design. It consists of two steps including LMC as well as GMC. While the local motion is specific to an individual moving object, the global motion, defined as a movement which has a common impact on all voxels, may originate from various aspects such as subject movement, maternal respiration, cardiovascular motion, and instrumental intervention (e.g., vibration).

The DOMC procedures for the hyperoxia studies of the placenta and fetal brain are shown in Fig. 3. The DOMC framework consists of two steps. The first step is to temporally split a whole volume sequence into independent phases and to define a phase-specific template which enables independent image registration at each phase. The phase-specific template can be chosen as either a single arbitrary volume inside the phase, or the mean volume averaged over the phase sequence in order to reduce the vulnerability of significant temporal variation of BOLD intensities, spin history artifact (induced by fetal motion), and data loss due to between-slice gap.³⁴ The second step is to restrict the FOV for image registration to the ROI and its surrounding adjacent regions in order to estimate the local motion of each moving object. This is in contrast to typical (or global) motion correction where the FOV is chosen to cover the full volume.³⁴ The ROI-restricted FOV can be obtained using the ROI masks which are usually defined from the phase-specific template.

Fig. 3

The structure of the DOMC. The DOMC consists of global and LMC steps. Each block corresponds to a volume consisting of voxel-wise BOLD signals at a specific time point.

All the above strategies can be mathematically formulated in a synthesized manner as follows. Let us consider an fMRI sequence $\mathbf{X}\equiv \{X_t;t=1,\dots ,N\}$ of length $N$ where $X_t$ is a three-dimensional volume and $X_t(\mathbf{q})$ denotes the intensity at voxel $\mathbf{q}$ in the volume $X_t$. Assume that $\mathbf{X}$ includes $K$ independently moving objects and can be separated into $\mathcal{P}$ independent phases according to a sequential block design. Then, the image registration is posed as the problem to estimate the combination of optimal global and local motion parameter set $[{\mathbf{m}}_G(t),{\mathbf{m}}_1(t),\cdots ,{\mathbf{m}}_K(t)]$ at time $t$ where ${\mathbf{m}}_G(t)$ and ${\mathbf{m}}_i(t)$ denote the global and local motion parameters corresponding to the $i$‘th ROI, respectively, and can be summarized as

2.4.

Outlier Rejection

2.4.1.

Volume outlier rejection

VOR is the process of eliminating significantly misaligned outlier volumes even after motion correction, which can be regarded as a complementary method for removing spurious data from misaligned volumes. There exist some automatic tools of evaluating volume outliers, including AFNI (3dTqual) and FSL (fsl_motion_outliers).^39–41 However, these tools do not adequately address small subject motion as well as the temporal variation of BOLD signals since volume outliers are detected based on spatial correlation between a motion-corrected volume and the template. In this section, we propose a method for VOR, using a probabilistic approach based on the temporal outliers of BOLD time series, which is optimized to be robust for temporal variation of signals.

Suppose that we have a series of volume outlier scores (VOS) as a measure of evaluating the accuracy of volume alignment; that is, $\rho =\{{\rho }_1,\cdots ,{\rho }_N\}$ where ${\rho }_t$ denotes the VOS at time $t$. The volume outliers can then be detected by thresholding the VOS values. In other words, the VOR is formulated as the problem of computing the missingness vector $\mathbf{g}=[g_1,\dots ,g_N]$ such that $g_t=0$ if ${\rho }_t<{\rho }_{th}$, otherwise $g_t=1$ where ${\rho }_{th}$ is a predefined threshold.

Let us consider the $k$‘th ROI ${\mathbf{F}}_k$ which has been wrongly aligned to a phase-specific template at time $t$; in other words, the motion parameter ${\widehat{\mathbf{m}}}_{k,t}$ has been wrongly estimated so that the $t$’th volume should be rejected as a volume outlier by setting the missingness as $g_t=1$. Let ${\mathbf{q}}_{k,0}\in {\mathbf{F}}_k$ be a voxel with fixed coordinates inside the $k$‘th ROI of the template (‘0’ denotes a template volume), and let ${\mathbf{q}}_{k,t}$ be the actual location at time $t$ corresponding to the voxel ${\mathbf{q}}_{k,0}$. The volume misalignment leads to the spatial mismatch between ${\mathbf{q}}_{k,0}$ and ${\mathbf{q}}_{k,t}$. For the sake of simplicity, the rotational components of volume misalignment are considered negligible. The geometric distance ${\mathbf{n}}_{k,t}$ between two positions is then equivalent to the measurement error of object motion; that is, ${\mathbf{n}}_{k,t}={\mathbf{q}}_{k,t}-{\mathbf{q}}_{k,0}={\mathbf{m}}_{k,t}-{\widehat{\mathbf{m}}}_{k,t}$. As a result, the actual area of the ROI at time $t$ is updated as ${\mathbf{F}}_{k,t}\equiv \{{\mathbf{q}}_{k,t}|{\mathbf{q}}_{k,t}={\mathbf{q}}_{k,0}+{\mathbf{n}}_{k,t},{\mathbf{q}}_{k,0}\in {\mathbf{F}}_k\}$, and consequently the background at time $t$ would be given as ${\mathbf{B}}_{k,t}\equiv \mathbf{U}-{\mathbf{F}}_{k,t}$ where $\mathbf{U}$ is a set of whole voxels. Finally, the volume misalignment leads to the spatial mismatch between ${\mathbf{F}}_k$ and ${\mathbf{F}}_{k,t}$ as illustrated in Fig. 4.

Fig. 4

Motion-induced overlapping of ROI between two subsequent images. Fetal motion produces a nonoverlapping part of ROI between two subsequent volumes. Large motion leads to the large nonoverlapping area, which means that the probability of temporal outliers in the nonoverlapping region increases. The MOP can be used as the VOS on the basis of the Bayes’ rule.

To represent it in terms of probability, let $P_q(B_{k,t})\equiv P({\mathbf{q}}_{k,0}\in {\mathbf{B}}_{k,t}|{\mathbf{n}}_{k,t})$ denote the probability of the voxel ${\mathbf{q}}_{k,0}$ not belonging to the updated ROI ${\mathbf{F}}_{k,t}$ in time $t$ due to the error motion ${\mathbf{n}}_{k,t}$. As long as the shape of ROI ${\mathbf{F}}_k$ is always convex, the mean of $P_q(B_{k,t})$ over the ROI voxels, denoted by ${\rho }_k(t)$, is a monotonically increasing function $f_R$ over $|{\mathbf{n}}_{k,t}|<|{\mathbf{n}}_{\max }|$; that is,

Eq. (2)

$${\rho }_k(t)\equiv \frac{1}{V_k}\sum _{{\mathbf{q}}_{k,0}}P({\mathbf{q}}_{k,0}\in {\mathbf{B}}_{k,t}|{\mathbf{n}}_{k,t})={\zeta }_R(|{\mathbf{n}}_{k,t}|),$$

where $|{\mathbf{n}}_{\max }|$ is the maximum geometric error, $V_k\equiv N({\mathbf{F}}_{k,0})$ is the number of all voxels in the $k$‘th predefined ROI, and ${\zeta }_q(·)$ is a monotonically increasing function over the degree of translation $|{\mathbf{n}}_{k,t}|$.^42,43 If we especially assume that the function ${\zeta }_R$ is approximately proportional to $|{\mathbf{n}}_{k,t}|$, it can be exploited as the VOS ${\rho }_t$ to evaluate the degree of translation. (Note: this is not exactly proportional in practice because of the dependence of $P_q(B_{k,t})$ on the object shape).

Let ${\rho }_k(t)$ be the MOP of the $k$‘th ROI at time $t$. The MOP can be computed by estimating $P_q(B_{k,t})$. According to the law of large numbers,⁴⁴ $P_q(B_{k,t})$ can be approximated as $V_k^{(b)}(t)/V_k$ where $V_k^{(b)}(t)=N({\mathbf{F}}_{k,0}\cap {\mathbf{B}}_{k,t})$ denotes the number of voxels in the $k$‘th ROI not belonging to the updated ROI due to motion in time $t$; therefore, ${\rho }_k(t)\approx V_k^{(b)}(t)/V_k$. While $V_k^{(b)}(t)$ is determined depending on ${\mathbf{n}}_{k,t}$, the geometric error is unknown. To estimate the MOP ${\rho }_k(t)$, the temporal variation of BOLD signals can be exploited as an alternative. As depicted in Fig. 4, the volume misalignment leads to temporal outliers in BOLD signals at some boundary voxels of the predefined ROI. The typical approach to detecting temporal outliers is to calculate the TOS based on the Euclidean distance of intensities from the median absolute deviation (MAD) of the BOLD signals.¹⁶ Let $S_q=[s_1,s_2,\cdots ,s_N]$ be a vector of TOSs corresponding to the BOLD signal ${\mathbf{x}}_q=[x_q(1),x_q(2),\dots ,x_q(N)]$ at the voxel ${\mathbf{q}}_{k,0}$. Then, the voxel $\mathbf{q}$ is regarded as an outlier at time $t$ if $s_t>s_o$ where the threshold $s_o$ is automatically determined as the voxel-wise IQR boundary, in other words, $s_o=1.5\times \mathrm{STD}(S_q)$.

Let us consider the probability of a boundary voxel having temporal outlier in its BOLD signal due to volume misalignment. Let $P_q(O_{k,t})\equiv P[\omega ({\mathbf{x}}_q,t)=O]$ be the outlier probability, that is, the probability of a voxel ${\mathbf{q}}_{k,0}$ belonging to the outlier class $O$ in time $t$ where $\omega (·)$ denotes the class of BOLD signal. Then, we have the following relationship according to the Bayes’ rule:⁴⁵

Eq. (3)

$$P_q(O_{k,t})=P_q(O_{k,t}|B_{k,t})P_q(B_{k,t})+P_q(O_{k,t}|F_{k,t})[1-P_q(B_{k,t})]\mathrm{.}$$

In a similar manner with $P_q(B_{k,t})$, $P_q(O_{k,t})$ can be approximated to $V_k^{(o)}(t)/V_k$ where $V_k^{(o)}(t)$ is the number of all outlier voxels in $\mathbf{U}$ in time $t$. From Eq. (3), we obtain

Eq. (4)

$${\rho }_k(t)\approx \frac{V_k^{(o)}(t)/V_k-P_q(O_{k,t}|F_{k,t})}{P_q(O_{k,t}|B_{k,t})-P_q(O_{k,t}|F_{k,t})}\mathrm{.}$$

As shown in Fig. 5(a), the posterior probabilities $P_q(O_{k,t}|B_{k,t})$ and $P_q(O_{k,t}|F_{k,t})$ can be computed by integrating the corresponding posterior probability density function (PDF) over the outlier region where the TOS $s_t$ is larger than the threshold $s_o$; for instance, $P_q(O_{k,t}|B_{k,t})$ is given as follows

Eq. (5)

$$P_q(O_{k,t}|B_{k,t})={\int }_{s_o}^{\infty }p(s|B_{k,t})\mathrm{d}s\mathrm{.}$$

The background posterior probability $P_q(O_{k,t}|B_{k,t})$ depends on how far apart two PDFs are from each other as well. Let us suppose that the BOLD intensities in the $k$‘th ROI $F_{k,t}$ have statistically significant differences from those of the surrounding background $B_{k,t}$ so that two density distributions $p(s_t|F_{k,t})$ and $p(s_t|B_{k,t})$ get further away from each other, as exemplified in Fig. 5(b). There exists a critical point $s_c$ of TOS as the boundary between two classes where all voxels with TOS $s_t<s_c$ are assigned to the $k$‘th ROI, otherwise to the background. As two PDFs become further apart from each other, $s_c$ becomes larger than TOS threshold $s_o$; that is, $s_o\le s_c$. As a result, most voxels in the background would have temporal outliers in time $t$ as illustrated in Fig. 5(a). Finally, the background posterior probability $P_q(O_{k,t}|B_{k,t})$ would have a high value close to 1.

Fig. 5

The class-conditional probability density functions of temporal outliers in the foreground and background. (a) A theoretical model is visualized where the PDFs $p(s|B_k)$ and $p(s|F_k)$ are assumed to be normally distributed over TOSs $s$. The decision of a voxel being foreground or background is made at the critical point $s_c$ where $p(s_c|F_k)P(F_k)=p(s_c|B_k)P(B_k)$. On the other hand, the decision of a voxel being a temporal outlier is made at the boundary $s_o$ corresponding to a predefined rate of the normal distribution. (b) The difference in probabilities of relative BOLD intensities, corresponding to outlier score $s$, is illustrated between the fetal brain and its surrounding background by using intensity histograms. Equation (5) is still valid although the intensity distribution of surrounding background is not perfectly normal.

Normally, ${\chi }_{ob}\equiv P_q(O_{k,t}|B_{k,t})$ and ${\chi }_{of}\equiv P_q(O_{k,t}|F_{k,t})$ can be determined as some high and small values which satisfy ${\chi }_{of}<\min V_k^{(o)}(t)/V_k$ and $\max V_k^{(o)}(t)/V_k<{\chi }_{ob}$ since ${\rho }_k(t)\in [\mathrm{0,1}]$ and we assume that ${\chi }_{ob}$ and ${\chi }_{of}$ are almost constant over time. Note that those probabilities do not need to be precisely estimated practically, since their inaccuracy can be compensated by adjusting the threshold ${\rho }_{th}$. Finally, the missingness vector $\mathbf{m}$, which describes whether each volume is trimmed or not, is produced by thresholding the MOP ${\rho }_k(t)$ in each volume.

2.5.

Data Imputation

The process of VOR produces missing data in ROI-averaged time series, which results in the discontinuity in BOLD signals, or irregular time intervals between two subsequent volumes. The problem can be overcome by imputing the missing data in a statistical manner.

Let us consider $K$ ROI-averaged signals of length $N$ which were averaged at each ROI; ${\tilde{\mathbf{y}}}_i={[{\tilde{y}}_i(1),\cdots ,{\tilde{y}}_i(N)]}^T$ denotes the mean BOLD signal of the $i$‘th ROI ($i=1,\cdots ,K$). Assuming that the signals are incomplete in that they have zero value at each missing data point, let ${\mathbf{M}}_i=\mathrm{diag}[g_{i,1},\cdots ,g_{i,N}]$ be the missingness matrix for the $i$‘th ROI where $g_{i,t}=1$ if observed at time $t$ and $g_{i,t}=0$ otherwise; it is simply the matrix extension of the missingness vector $\mathbf{g}$. Let ${\mathbf{y}}_i$ denote the complete data without missing data corresponding to ${\tilde{\mathbf{y}}}_i$, and consider the regression model ${\mathbf{y}}_i={\mathbf{f}}_i+{ϵ}_i$ where $f_i(t_0)=E(y_i|t=t_0)$ is a finite continuous function and ${ϵ}_i(t)$ is an independent and identically distributed random variable with zero mean and variance ${\sigma }_i^2$, that is, ${ϵ}_i\sim N(0,{\sigma }_i^2)$. Our goal is to estimate the complete data ${\mathbf{y}}_i$ from the incomplete data ${\tilde{\mathbf{y}}}_i={\mathbf{M}}_i{\mathbf{y}}_i$ by imputing missing data as follows

Fig. 6

An example of data imputation based on local polynomial regression. The derivatives at a time point are obtained by applying a kernel to its neighboring data inside the predefined bandwidth of finite size. All missing data points are estimated based on the regressed time series.

2.6.

Performance Evaluation

The performance of motion correction was evaluated using the voxel-wise mean absolute residual variation (m-ARV) score as well as the ratio of temporal outliers in the resulting BOLD signals, and was compared with the FMRIB’s Linear Image Registration Tool (FLIRT)³⁴ which has been one of the conventional motion correction tools used in fMRI studies. The m-ARV score $r_q$ for the voxel $\mathbf{q}$ is defined as the absolute mean value of residual intensity variations over time; that is, $r_q=N^{-1}\sum _{t=1}^N|X_t(\mathbf{q})-\overline{X}(\mathbf{q})|$ where $\overline{X}$ is the $t$-mean volume over a phase.³⁴

While the linear registration method has been used as the default in each phase sequence for each ROI through the proposed DOMC pipeline, we also tested the effects of advanced image registration methods specialized for fetal motion correction. The test methods were made up of a possible combination of slice-to-volume registration and nonlinear image registration which can be regarded as two key features of the existing methods for fetal motion correction.^13,48 One of the simplest methods for slice-to-volume registration was implemented according to the slice timing correction algorithm, proposed in Ref. 48, which decomposes each volume into two subvolumes of odd and even slices, and aligns two subvolumes to remove the spatial mismatch between slices due to interleaved data acquisition. On the other hand, the FMRIB tool for small-displacement non-linear registration (FNIRT) was utilized for nonlinear image registration.⁴⁹

To determine the impact of either fetal motion or acquisition plane on the quality of fMRI data, we analyzed the association between estimated motion parameters and the ratio of rejected volume outliers (RVO), defined to be the fraction of rejected volumes to total volumes. The performance of the proposed LPS-based data imputation was evaluated through the random attack experiment where some data are randomly removed in a time series, and quantified using the data imputation score (DIS) which is defined as the correlation coefficient between ground truth and recovered time series. To compare the statistics of estimated motion parameters either between different temporal phases or between healthy and CHD fetuses, the absolute distance of translation (ADT) was defined as the norm of a translation vector.

3.1.

Design-Optimized Motion Correction

The quality of each EPI sequence was first visually inspected using the FSL viewer before preprocessing. Four subjects were found to have several erroneous slices due to either Nyquist ghosting or slice shifting effects. These erroneous slices were excluded from image registration.

Figure 7 shows examples of fetal fMR images after applying either the proposed DOMC method or FLIRT. The top row images depict examples where the motion of the fetal brain was not well aligned to the actual location when FLIRT was applied for motion correction. On the other hand, the image registration for the fetal brain was improved using the proposed DOMC method as shown in the bottom row, where the registered ROIs coincided with the actual corresponding regions.

Fig. 7

Comparison of DOMC with FLIRT in the fetal brain. The fMRI image registered by DOMC (b) was spatially well matched with the predefined ROI area (highlighted by red) while the image registered by FLIRT (a) was not as indicated by arrows.

Figure 8 illustrates the effect of the proposed DOMC method and FLIRT on a ROI-averaged BOLD time series. The averaged time series were corrupted by many temporal outliers caused by volume misalignment after applying the standard preprocessing pipeline based on FLIRT. Once the proposed DOMC method was applied, the number of temporal outliers in the averaged BOLD time series was reduced.

Fig. 8

The effect of the proposed DOMC method on an ROI-averaged BOLD time series of the fetal brain, compared to the standard motion correction approach based on FLIRT. Some temporal outliers induced by fetal motion (indicated by arrows) were prominently reduced when the proposed DOMC method was applied.

Figure 9 summarizes the quantitative comparison of performance in motion correction between the proposed DOMC method and FLIRT. The ratios of temporal outliers in the resulting BOLD signals after motion correction were significantly reduced using the proposed method compared to FLIRT in both the fetal brain ($1.18\pm 1.36$ vs. $4.55\pm 3.03$, $p<0.0001$ using the paired one tail t-test) and placenta ($0.16\pm 0.37$ vs. $0.78\pm 0.83$, $p=0.009$).

Fig. 9

The ratios of temporal outliers in the mean BOLD signals of the brain and placenta comparing the DOMC and FLIRT. The asterisk symbol indicates that the group difference in the ratio of temporal outliers is statistically significant with $p$ value $<0.05$.

Figure 10 shows that the mean m-ARV score exhibited the tendency of decreasing progressively over the sequential steps of GMC and LMC in both the fetal brain and placenta, which indicates that the similarity between volumes was improved. Figure 10 also shows the quantitative comparison of motion correction performance between five different image registration methods applied to the DOMC pipeline. The nonlinear registration approach consistently exhibited more enhanced performance of motion correction, with lower m-ARV scores ($28.6\pm 5.5$ vs $38.7\pm 4.7$ for the fetal brain, $24.5\pm 1.1$ vs $44.4\pm 3.3$ for the placenta), than the linear registration in both the fetal brain and placenta. However, its computation speed was approximately 2.62 times slower (729 vs 278 minutes on average). On the other hand, the effect of slice-to-volume registration on image quality was conspicuous only in the fetal brain when it was followed by nonlinear motion correction (N vs NS: $25.4\pm 3.8$ vs $31.8\pm 5.5$ for m-ARV), while its effects were negligible in the placenta. Regardless of image registration methods, the motion correction was always more effective in the baseline compared to hyperoxia or return to baseline.

Fig. 10

The effects of nonlinear image registration and slice-to-volume registration applied to the common DOMC pipeline on the motion correction performance. Lower m-ARV scores correspond to better image quality after motion correction. GMC denotes applying the GMC only while LMC denotes applying both the global and LMCs. For the latter case, we compared four possible variations: linear motion correction (L) which was set as the default of the DOMC pipeline, linear motion correction along with slice-to-volume registration (LS), nonlinear motion correction (N), and nonlinear motion correction along with slice-to-volume registration (NS).

3.2.

Outlier Rejection and Data Imputation

Figure 9 quantitatively demonstrates that some outlier volumes may still exist due to the failure in estimating the motion parameters. Figure 11 illustrates a detailed example which shows how the failure in motion correction leads to a volume outlier. The top row shows two consecutive images with a large rotation of the fetal brain resulting in a considerable difference in the shape and location of the brain projected to a specific plane. This resulted in a number of temporal outliers especially on the boundary voxels of the fetal brain. The TOS was obtained in each voxel, and the sum of TOS over all ROI voxels was computed and converted into the MOP as a VOS. As shown in the bottom row graph, the VOS at the $(t+1)$’th volume was greater than the threshold which was set as 0.3, and therefore the volume was rejected.

Fig. 11

An example of VOR. The high motion between two images causes temporal outliers in the boundaries of the fetal brain and finally leads to the substantial increase in the VOS.

Table 2 shows the dependence of RVO according to types of acquisition plane. The RVO was consistently higher in axial view compared to the coronal view, regardless of subject types. Figure 12 shows an example of detected outlier voxels in the placenta before averaging voxel-wise BOLD signals, demonstrating that the (red-highlighted) voxels corrupted by fetal body or nonplacental tissues were automatically detected as outliers.

Table 2

The ratios of rejected volume outliers (RVO) in coronal and axial planes.

Fig. 12

An example of voxel outliers in the placenta. The original ROIs were shaded with yellow color while the voxel outliers were shaded with red color.

Figure 13 shows the comparison of the proposed LPS-based data imputation method (based on local polynomial smoothing) with other traditional methods based on global curve fitting which was applied in Ref. 31. In the LPS-based data imputation, the DIS was maximized over the whole range of missing volumes in both the fetal brain and placenta when the bandwidth of local smoothing was set to be 20. In the fetal brain, the DIS gradually decreased as the percentage of missing data increased. On the other hand, the DIS in the placenta was maintained over 0.95 until 30% of the data were missing while it sharply declined when more data than 30% were missing. The traditional method based on the Fourier series model exhibited similar performance to the LPS-based method; however, the DIS was slightly improved, particularly with 15% to 30% missing data, using LPS-based data imputation compared to the Fourier series model.

Fig. 13

The comparison of performance among data imputation techniques both in (a) fetal brain and (b) placenta. Some time points were randomly eliminated according to a controlled percentage of missing data, and were recovered using a diversity of data imputation techniques. Then, the correlation coefficient between the original ROI time series (as a ground truth) and the recovered time series was computed to evaluate the performance of data imputation. Both ‘Poly1’ and ‘Fourier3’ denote the linear polynomial model and the Fourier series model with three degrees, respectively, and $W$ denotes the bandwidth of local polynomial smoother with the polynomials with degree 5.

3.3.

Clinical Application

Table 3 summarizes the statistics of estimated motion parameters in different phases for both healthy and CHD fetuses. The mean ADTs were smaller in GMC than in LMC regardless of phases and subject type (i.e., controls and CHD cases). On the other hand, the mean ADTs were higher in CHD fetuses compared to healthy control fetuses for both global and LMC; however, the differences in mean ADTs were not statistically significant. Specifically, the mean ADTs of the fetal brain were 20.7, 51.1, 14.7 times those of the placenta at baseline, hyperoxia, and return to baseline, respectively. Interestingly, the standard deviation of ADTs for the fetal brain of CHD cases was larger than that for healthy controls in all phases.

Table 3

The statistics of absolute translations estimated by global and LMC. P1, P2, and P3 correspond to baseline, hyperoxia, and return to baseline, respectively.

The phase-dependent mean ADTs of global motion were commonly increased during both hyperoxia and return to baseline compared to baseline before oxygen supply. In the LMC, the phase-dependent changes in the mean ADTs of local motion were considerably different between healthy controls and CHD fetuses. In the healthy controls, the mean ADTs of the fetal brain and placenta were maintained or slightly increased during both phases of hyperoxia and return to baseline. Conversely, the mean ADTs of brain motion in CHD fetuses were decreased after hyperoxia.

The group differences in RVO between healthy and CHD fetuses are shown in Table 2; the RVO was reduced during hyperoxia in healthy controls, while it was increased in fetuses with CHD. The RVO was normally higher in CHD fetuses compared to healthy control fetuses.

4.1.

Methodological Issues

This paper proposed a preprocessing pipeline to overcome the practical limitation of fetal motion correction and eliminating significant noise artifacts that are unique challenges in stimulus-based fMRI of the living fetus. The novelties of this paper can be summarized as follows. First, motion correction is optimized to the experimental design. Second, the misaligned volumes are automatically detected using a probabilistic approach optimal to fMRI data, and missing data resulted from VOR are recovered through data imputation to allow for advanced and reliable BOLD time-series analysis retrospectively.

Through the analyses of motion parameters estimated by the proposed DOMC method, we also found that the pattern of fetal motion changes over different stimulus phases.⁵⁰ Fetal motion tended to be dependent on the task design where the degree of fetal motion changed over the study phases. In addition, the degree of absolute translation was significantly smaller in the placenta compared to the fetal brain by visual inspection of EPI sequences. These results support the need for DOMC, based on the assumption of heterogeneous patterns of movement over different phases and ROIs.

The different pattern of movement over phases does not matter in resting state fMRI where no stimulus is assigned. Since an EPI sequence acquired at rest does not need to be temporally split into several phases, the DOMC method becomes roughly identical to the typical image registration which has been applied in other fMRI studies on the fetal brain.^2,26,27 On the other hand, if the EPI sequence is acquired at rest but contains multiple organs moving independently, it still needs to be split into different ROIs so that the DOMC method can be applied. As such, the proposed DOMC method may also be beneficial for resting state fMRI of the moving fetus.

The performance of the proposed DOMC method was compared with FLIRT. We showed (Fig. 9) that the number of outliers was significantly reduced in the resulting averaged BOLD signals after DOMC. Although the reduction of temporal outliers in averaged BOLD signals is not a direct metric for performance evaluation, it might have been indirectly attributed to the improved performance of motion correction, assuming that the volume outliers are dominated by high fetal motion.

We also compared four advanced variations of image registration in the proposed DOMC pipeline with linear image registration to determine the effect of nonlinear image registration and slice-to-volume registration. Our results demonstrate that nonlinear image registration leads to better performance in motion correction than linear registration in both the fetal brain and placenta. This is not surprising, particularly in the placenta which is nonrigid but deformable as shown in Fig. 12. Nevertheless, we used linear image registration in the present study to minimize the computational complexity. On the other hand, the slice-to-volume registration did not always lead to a substantial improvement in motion correction, with the exception of nonlinear registration when applied to the fetal brain. This negligible effect might be attributed to the simplicity of the applied slice-to-volume motion correction algorithm in which two sub-datasets of odd and even slices are simply realigned. More advanced methods, as proposed in Ref. 13, need to be developed to improve the accuracy of slice-to-volume registration for both the fetal brain and placenta. Taken together, these results demonstrate that the existing methods for motion correction could be effectively incorporated into the proposed DOMC pipeline to improve the performance of motion correction although their high computational complexity needs to be enhanced for practical use.

The performance of LPS-based data imputation was quantitatively evaluated through the random attack experiments. We showed (Fig. 13) that the performance is highly dependent to the bandwidth of local smoothing. While the bandwidth was empirically determined in our study, it would be valuable to develop an algorithm to determine the optimal bandwidth automatically according to the statistical properties of individual time series to maximize the performance of data imputation compared to other traditional methods. Interestingly, the data imputation was more successful in the placenta than in the fetal brain as long as the missing rate is kept low ($<30\%$).

In this study, we compared coronal and axial acquisitions to examine the impact of the acquisition plan on the performance of motion correction. The number of rejected volume outliers in the fetal brain was higher with axial acquisitions compared to coronal acquisitions. These data demonstrate that the motion correction of the fetal brain is less effective in the axial plane than in the coronal plane. This effect may be attributed to the interleaved data acquisition given that the interleaved scan resulted in spatial mismatch between slices due to fetal motion and finally led to the inaccuracy in estimating motion parameters relevant to the between-slice movement. In the case of axial acquisitions, the performance of motion correction may be deteriorated since the EPI sequence includes more between-slice movement of the fetal brain while it appropriately captures the in-plane rotation of the brain. On the other hand, in the coronal plane, motion correction was more successful in estimating the translational movement of fetuses and less accurate in detecting the rotation of the fetal brain which is normally related to the between-slice motion in EPI volume. Our data suggest that the coronal acquisition may be more beneficial than the axial acquisition resulting in improved performance of motion correction of the fetal brain.

Conversely, the acquisition plan had a different impact on motion correction of the placenta. The placenta demonstrated periodic motion attributed to maternal respiration. As observed empirically in our study, the anterior-posterior components of respiration-related placental motion artifact could be mitigated using data acquired in the axial plane compared to the coronal plane, which may be attributed to the anthropometric features of the placenta which is a flat organ and is typically attached to the lateral wall of the uterus. However, this conflicts with our conclusion that the coronal plane is better for the fetal brain. Since placental motion was not as pronounced compared to the fetal brain, we chose the coronal plane as a trade-off to correct more serious motion of the fetal brain. Indeed, the changes in the number of rejected volume outliers were fewer in the placenta.

Although there are a number of strengths to our preprocessing pipeline, the limitations also deserve mention. The first limitation is that a unique set of parameters in computing of VOS for VOR were uniformly applied to all subjects; however, they may be different over subjects or over ROIs. The VOS is determined depending on background and foreground posterior probabilities [as described in Eq. (4)] which may be different according to either ROI or subject types. For example, the temporal outliers tend to frequently occur in the inner local part of an ROI with spatially inhomogeneous intensities. In other words, an ROI voxel is more likely to have a temporal outlier in the case that the ROI is highly textured, which leads to the increase in the foreground posterior probability. Moreover, those posterior probabilities have been approximately established under such an inherent assumption that each ROI has significantly different BOLD intensities compared to the surrounding background. While this assumption is broadly appropriate for most ROIs including the placenta, fetal liver, and heart, it does not always hold true for the fetal brain which tends to have fewer intensity differences compared to the surrounding amniotic fluid and uterine tissues.¹³ Therefore, an automatic algorithm for deciding the threshold for VOS adaptively in each subject or each ROI would be instrumental to enhance the accuracy in computing the VOS and ultimately improve the performance of VOR. Another limitation is that our proposed preprocessing pipeline does not directly address physiological artifacts such as cardiac pulsation and respiratory motion which may have considerable effects on fetal motion correction and missing data imputation. Although the physiological artifacts could be compensated in part by the outlier rejection process, an advanced method of eliminating such physiological noises found in fetal MRI will need to be developed.

4.2.

Clinical Application to Congenital Heart Disease

We showed that the proposed preprocessing pipeline can be effectively employed to characterize fetal motion in healthy controls and CHD fetuses. Our preliminary data suggest that the degree of fetal motion tends to increase during hyperoxia in CHD fetuses (but not significantly). In addition, the motion of the fetal brain in CHD cases showed higher variance during hyperoxia compare to controls. These observations suggest that the CHD fetus may be more responsive to maternal hyperoxia. However, these pilot data need to be validated on a larger cohort of healthy and high-risk CHD fetuses. Investigating the mechanisms underlying the differences observed between these two groups was beyond the scope of this study but currently underway.

4.3.

Conclusion

In conclusion, the proposed preprocessing pipeline including DOMC and outlier rejection can be effectively used to eliminate noise artifacts induced by fetal and placental movement in stimulus-based fMRI. Although the pipeline needs to be improved through further studies to ensure accurate and reliable signal analyses, this work makes an important technical advance for robust preprocessing in stimulus-based functional neuroimaging studies of the fetal brain and placenta, and lays the foundation not only for noninvasive functional assessment of fetal brain-placenta unit, but also for enabling early detection of impaired fetal brain-placenta circulation in future.