Home | About | Journals | Submit | Contact Us | Français |

**|**J Med Imaging (Bellingham)**|**v.3(2); 2016 April**|**PMC4821892

Formats

Article sections

Authors

Related links

J Med Imaging (Bellingham). 2016 April; 3(2): 026001.

Published online 2016 April 5. doi: 10.1117/1.JMI.3.2.026001

PMCID: PMC4821892

Wonsang You,^{a}^{} Iordanis E. Evangelou,^{a,}^{b}^{} Zungho Zun,^{a}^{} Nickie Andescavage,^{c,}^{d,}^{e}^{} and Catherine Limperopoulos^{a,}^{b,}^{c,}^{e,}^{*}^{}

Received 2015 August 31; Accepted 2016 March 3.

Copyright © The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Fetal motion manifests as signal degradation and image artifact in the acquired time series of blood oxygen level dependent (BOLD) functional magnetic resonance imaging (fMRI) studies. We present a robust preprocessing pipeline to specifically address fetal and placental motion-induced artifacts in stimulus-based fMRI with slowly cycled block design in the living fetus. In the proposed pipeline, motion correction is optimized to the experimental paradigm, and it is performed separately in each phase as well as in each region of interest (ROI), recognizing that each phase and organ experiences different types of motion. To obtain the averaged BOLD signals for each ROI, both misaligned volumes and noisy voxels are automatically detected and excluded, and the missing data are then imputed by statistical estimation based on local polynomial smoothing. Our experimental results demonstrate that the proposed pipeline was effective in mitigating the motion-induced artifacts in stimulus-based fMRI data of the fetal brain and placenta.

Advances in *in vivo* functional magnetic resonance imaging (fMRI) of the fetus are improving our understanding of fetal brain development in healthy and high-risk pregnancies.^{1}^{}^{–}^{3} Similar to traditional fMRI studies,^{4} fetal fMRI studies are performed by measuring either stimulus-based activity or spontaneous activity at rest. Recently, changes in the oxygenation of the fetal brain and placenta during maternal hyperoxia have been explored using blood oxygenation level dependent (BOLD) fMRI to study the compromised fetus *in utero*.^{5}^{,}^{6}

However, such hyperoxia studies are technically challenged by degradation of the BOLD signal attributed primarily to fetal and maternal movement. The fetus often moves continuously during the MR study^{7} while the adult has the tendency to exhibit slow drift or intermittent spike-like head motion.^{8}^{}^{–}^{10} Consequently, fetal motion has a significant impact on BOLD intensities given that the signal depends on both the spatial position^{8} and the spin excitation history of magnetization.^{11}^{}^{–}^{13} This motion becomes more problematic in stimulus-based fMRI than in stimulus-free fMRI, given that the fetus tends to move more in response to the stimulus, which in turn causes serious spatial misalignment between volumes at different time points.^{14}^{,}^{15} Moreover, each fMRI volume includes both the fetal brain and placenta, which exhibit different types of motion. The fetal brain behaves as a rigid body with high range of motion, while the placenta exhibits a passive, low range of motion—but because of its flexible underlying structure, it often undergoes deformations due to fetal and maternal movement. The placenta may also exhibit nonisotropic motion due to maternal breathing, local contraction and relaxation, and movement of tissues surrounding the uterine cavity.

In addition to motion-induced artifacts, a number of other artifacts of physiological origin perturb the BOLD fMRI signals. While breathing and cardiac pulsation are considered major non-neuronal sources of physiological noise in the adult brain, the BOLD signals in fetal fMRI are also strongly influenced by other physiological artifacts originating from amniotic fluid flow and maternal bowel movement in the abdominal cavity. Air or iron in the maternal intestines may also lead to significant signal loss in regions of interest (ROI), the so-called bubble artifact.^{13}

These artifacts observed in stimulus-based fMRI of the fetus limit the success of traditional approaches for motion correction that are based on rigid-body image registration (implemented in AFNI,^{16} AIR,^{17}^{,}^{18} FSL,^{19} and SPM^{20}^{,}^{21}). This is in large part due to the fact that these approaches make assumptions such that the MR image data contain either the brain or a single rigid-body object, and that the background intensity variation is negligible.

To date, technical advances in motion correction for MR images of the moving fetus^{22}^{}^{}^{–}^{25} have focused on anatomical MR images of the fetal brain. Very few studies have addressed fetal motion correction in fMRI studies.^{1}^{,}^{2}^{,}^{13}^{,}^{26}^{,}^{27} Fulford et al.^{28}^{,}^{29} eliminated maternal signals outside the fetal skull to avoid the independent movement between maternal and fetal parts, and then applied conventional image registration using visual inspection to stimulus-based fMRI of the fetal brain. Ferrazzi et al.^{13} suggested an MR-based preprocessing pipeline for resting-state fMRI of the fetus, which consisted of bias field correction, slice-to-volume registration, and spin history correction. More recently, Seshamani et al.^{12} suggested a method based on a bias field model to correct intensity inhomogeneity caused by fetal motion. Despite these technical advances, motion correction methods for fMRI studies of multiple organs of the fetus remain sparse. Less than a handful of studies on hyperoxia fMRI data of the placenta and fetal organs have addressed the motion artifact problem;^{5}^{,}^{6} however, they attempted to segment ROIs manually in each volume to circumvent the practical limitation of typical motion correction in such a complex environment.

While available fetal fMRI studies have commonly adopted traditional image registration with minor variations for fetal motion correction, they have shown that the typical approach to motion correction was not always successful due to instantaneous high motion, which led to a considerable number of volumes remaining significantly misaligned after motion correction was applied.^{2}^{,}^{26}^{,}^{27} Those misaligned volumes were identified either manually through visual inspection or automatically by thresholding the mean voxel displacement, and were indispensably excluded as volume outliers so that only the “surviving” volumes could be exploited for postprocessing analyses. One important issue in volume outlier rejection (VOR) is that the typical criteria enforced for identifying a volume as an outlier are not ideally suited for fetal fMRI data. For example, one study attempted to remove such contaminated volumes manually, which reliability might be limited by subjective visual inspection used to assess the quality of volumes.^{1}^{,}^{30} The other issue is that VOR ultimately leads to irregularly distributed missing data points in voxel-wise BOLD time series, which markedly restricts the scope of available postprocessing data analyses.

Therefore, the objective of this paper was to develop a robust preprocessing pipeline dedicated to stimulus-based fMRI containing more than one independently moving ROI such as the placenta and fetal brain. We have previously introduced our initial strategies using a robust motion correction method for fMRI data of the fetal brain and placenta acquired with a single block paradigm of maternal hyperoxia.^{31} This paper presents the underlying mathematical foundations in which these preliminary strategies are extensively applied to more extensive stimulus block designs with slow cycle of phase transition (i.e., the phase length $>30\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$) as well as multiple ROIs of the moving fetus. In addition, we report an initial clinical application of our proposed method, by analyzing the group differences in fetal motion between healthy fetuses and fetuses diagnosed with complex congenital heart disease (CHD).

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\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{2}$ or $440\times 440\text{\hspace{0.17em}\hspace{0.17em}}{\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.

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).

The hyperoxia task paradigm consisted of three phases: a 2-min normoxia (baseline, 21% oxygen), followed by 100% oxygen ($15\text{\hspace{0.17em}}\mathrm{L}/\mathrm{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.

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.^{12}

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).^{34} The registration parameters were as follows: normalized cross-correlation ratio cost function, 256 histogram bins, trilinear interpolation, and search angle of $-20\text{\hspace{0.17em}\hspace{0.17em}}\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.^{35} 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.

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.^{34} 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.^{34} The ROI-restricted FOV can be obtained using the ROI masks which are usually defined from the phase-specific template.

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

$${({\widehat{\mathbf{m}}}_{G},{\widehat{\mathbf{m}}}_{1},\dots ,{\widehat{\mathbf{m}}}_{K})}_{t}=\underset{{\mathbf{m}}_{G},{\mathbf{m}}_{1},\dots ,{\mathbf{m}}_{K}}{\mathrm{argmin}}\text{\hspace{0.17em}}\sum _{k\in \mathrm{\Psi}}\sum _{\mathbf{q}\in {{\mathbf{R}}^{\prime}}_{k,i}}\xi ({X}_{t}(\mathbf{q}),{T}_{k,i}\{L[G(\mathbf{q};{\mathbf{m}}_{G});{\mathbf{m}}_{k}]\}),$$

(1)

where $k\in \mathrm{\Psi}=[1,\dots ,K]$, ${T}_{k,i}$ denotes the template of the $k$’th ROI in the $i$’th phase, $G$ and $L$ denote the geometric transformations of global and local motion, respectively. In the DOMC framework, the global and local motions can be separately estimated in a hierarchical manner. In the $i$‘th phase sequence, the motion parameter ${\mathbf{m}}_{k}(t)$ of the $k$’th ROI at time $t$ can be estimated by localizing the cost function $\xi $ to the voxels inside the dilated ROI mask ${{\mathbf{R}}^{\prime}}_{k,i}$.^{38} This in turn prevents the cost function from being stuck in local minima by excluding other moving objects and the dynamically changing background from the search field.^{34}

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.

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 **...**

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}}_{\mathrm{max}}|$; that is,

$${\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}|),$$

(2)

where $|{\mathbf{n}}_{\mathrm{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}(\xb7)$ 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,^{44}
${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.^{16} 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 (\xb7)$ denotes the class of BOLD signal. Then, we have the following relationship according to the Bayes’ rule:^{45}

$${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{.}$$

(3)

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

$${\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{.}$$

(4)

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

$${P}_{q}({O}_{k,t}|{B}_{k,t})={\int}_{{s}_{o}}^{\infty}p(s|{B}_{k,t})\mathrm{d}s\mathrm{.}$$

(5)

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.

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 **...**

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}<\mathrm{min}\text{\hspace{0.17em}}{V}_{k}^{(o)}(t)/{V}_{k}$ and $\mathrm{max}\text{\hspace{0.17em}}{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.

The BOLD intensities may be affected by various artifacts such as imprecise motion correction, physiological noises, and unexpected MR artifacts. For example, the limbs of the fetus may instantaneously deform the placenta due to its intense motion. Since it changes the local shape of the placenta, it is not easily captured by either motion correction or outlier rejection. A second example would be a gas bubble artifact which can be frequently introduced due to the maternal intestines surrounding the womb, which may result in significant signal loss in the local area of either the fetal brain or placenta. Such voxel outliers can be detected and rejected using the generalized extreme studentized deviate (ESD) procedure^{46} which is appropriate, compared to other typical methods such as Grubb’s test and Dixon’s test. Refer to Appendix A for its theoretical descriptions.

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}+{\u03f5}_{i}$ where ${f}_{i}({t}_{0})=E({y}_{i}|t={t}_{0})$ is a finite continuous function and ${\u03f5}_{i}(t)$ is an independent and identically distributed random variable with zero mean and variance ${\sigma}_{i}^{2}$, that is, ${\u03f5}_{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

$${\widehat{\mathbf{y}}}_{i}={\mathbf{M}}_{i}{\tilde{\mathbf{y}}}_{i}+({\mathbf{1}}^{(N)}-{\mathbf{M}}_{i})({\widehat{\mathbf{f}}}_{i}+{\widehat{\u03f5}}_{i}),$$

(6)

where ${\mathbf{1}}^{(N)}$ is the $N\times N$ matrix whose all entries are 1, and the error ${\widehat{\u03f5}}_{i}$ satisfies ${\widehat{\u03f5}}_{i}\sim N(0,\mathrm{var}[{\mathbf{v}}_{i}])$ with ${\mathbf{v}}_{i}\equiv \{x|x\in {\mathbf{M}}_{i}({\tilde{\mathbf{y}}}_{i}-{\widehat{\mathbf{f}}}_{i}),x\ne 0\}$. It indicates that the missing data can be obtained from the estimated regression function ${\mathbf{f}}_{i}$. A nonparametric estimator called the local polynomial smoother (LPS) can be effectively used to estimate the regression function ${\mathbf{f}}_{i}$ for diverse types of complex time series.^{47} In the LPS-based regression, the low degree polynomial is fitted, through weighted least squares, not to the whole time series but to temporal segments localized by a kernel (or moving window) with specific bandwidth (or smoothing parameter) as depicted in Fig. 6. Refer to Appendix B for mathematical descriptions.

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)^{34} 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.^{34}

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.^{49}

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.

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.

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.

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 **...**

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$).

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 $<$ **...**

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.

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.

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.

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.

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.

The statistics of absolute translations estimated by global and LMC. ${P}_{1}$, ${P}_{2}$, and ${P}_{3}$ 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.

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.^{50} 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.^{13} 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.

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.

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.

This work was supported by the Canadian Institutes of Health Research (MOP-81116) and the National Institutes of Health (RO1HL116585, IDDRC, grant P30HD40677). We thank Dr. Sonia Dahdouh for advice on theoretical description, and Dr. Josepheen Cruz for her input on time series analysis. We are very grateful to our research coordinators and the families that participated in our study.

•

**Wonsang You** is a research associate in the Developing Brain Research Laboratory at Children’s National Health Systems. He received his MS in engineering from the Korea Advanced Institute of Science and Technology in 2008, and his PhD in information technology from Otto-von-Guericke University Magdeburg in 2013. From 2009 to 2012, he worked for resting state fMRI in Leibniz Institute for Neurobiology. His research interests include in vivo neuroimaging data analyses of the fetal brain and placenta, and functional brain development of the fetus.

•

**Iordanis E. Evangelou** received his DPhil in medical imaging from the University of Oxford and his MS in management systems. He was a research scientist (2005 to 2011) at the National Institutes of Health, developing brain MR image acquisition and postprocessing. From 2011 to 2015, he was a MR scientist at Children’s National Medical Center and assistant professor of radiology and pediatrics at George Washington University, where he developed fetal brain MR spectroscopy methods and pediatric RF coils for brain and body imaging.

•

**Zungho Zun** is on the research faculty in the Developing Brain Research Laboratory at Children’s National Health System, and assistant professor of pediatrics at George Washington University. His main expertise is development of MRI pulse sequence and image reconstruction algorithms. His research interests lie in developing methodology of assessing hemodynamics in humans with the ultimate goal of clinical translation. Currently, his research mainly focuses on the development and validation of measuring perfusion in the placenta and fetal/neonatal brain.

•

**Nickie Andescavage** received her medical degree from Columbia University. Currently she is an attending neonatologist at Children’s National Health Systems, and an assistant professor of pediatrics at the George Washington University School of Medicine in Washington, DC. Her research interests include fetal brain development in healthy and complex pregnancies, as well as the early detection of brain injury in preterm and term infants through the use of advanced magnetic resonance imaging.

•

**Catherine Limperopoulos** is the director and principal investigator of the Developing Brain Research Laboratory at Children’s National Health Systems and associate professor of pediatrics at George Washington University School of Medicine and Health Sciences. Her research activities focus on studying the causes and consequences of early life brain injury in high-risk fetal, preterm, and full-term infant populations. Her long-term goals are to develop reliable biomarkers of brain injury that will guide medical and rehabilitation interventions aimed at circumventing injury.

Let us consider the $k$’th ROI ${\mathbf{F}}_{k}$ containing $Q$ voxels. Let ${\mathbf{F}}_{k}^{(s)}=\{{\mathbf{q}}_{1},{\mathbf{q}}_{2},\cdots ,{\mathbf{q}}_{Q}\}$ be the set of ordered voxels such that ${R}_{i}>{R}_{i+1}$ where ${R}_{i}$ is the absolute $z$-score for the $i$’th ordered voxel ${\mathbf{q}}_{i}$ with intensity $\mathbf{x}({\mathbf{q}}_{i})$; i.e., ${R}_{i}=|\mathbf{x}({\mathbf{q}}_{i})-\mu |/\sigma $ for the mean $\mu $ and standard deviation $\sigma $ of intensities over the voxel set ${\mathbf{F}}_{k}^{(s)}$. Then, a voxel ${\mathbf{q}}_{i}$, whose index $i$ is less than a predefined maximum, is considered as an outlier if its absolute $z$-score ${R}_{i}$ is larger than its corresponding critical value ${\lambda}_{i}$; i.e., ${R}_{i}>{\lambda}_{i}$ where

$${\lambda}_{i}=\frac{{t}_{(b,Q-i-1)}(Q-i)}{{(Q-i-1+{t}_{(b,Q-i-1)}^{2})}^{1/2}{(Q-i+1)}^{1/2}},$$

(7)

${t}_{(b,t)}$ is the $b$’th percentile of the Student’s $t$-distribution with $t$ DOF, and $b=1-\alpha /(Q-i+1)$ for the significance level $\alpha $ (which was set to be 0.05 in this study).

Assuming that the $(p+1)$’th derivatives of ${\mathbf{f}}_{i}$ at any $t$ exist, the regression function ${\mathbf{f}}_{i}$ can be estimated by finding the polynomial parameter vector ${\alpha}_{i}(t)={[{\alpha}_{i,0}(t),\cdots ,{\alpha}_{i,p}(t)]}^{T}$ where ${\alpha}_{i,j}(t)={f}^{(j)}(t)/(j!)$; in other words,

$${\widehat{f}}_{i}(t)={\widehat{\alpha}}_{i,0}(t)={\mathbf{e}}_{1}^{T}{\widehat{\alpha}}_{i}(t),$$

(8)

where ${\mathbf{e}}_{1}$ is a binary vector of length $p+1$ having 1 in the first entry and 0 for the rest. The parameter vector ${\alpha}_{i}(t)$ is estimated by minimizing the kernel-weighted local least squares as follows:

$${\widehat{\alpha}}_{i}(t|\mathbf{W},{\mathbf{M}}_{i})=\underset{{\alpha}_{i}(t)}{\mathrm{argmin}}\text{\hspace{0.17em}}{[{\tilde{\mathbf{y}}}_{i}-\mathbf{P}(t){\alpha}_{i}(t)]}^{T}\mathbf{W}(t){\mathbf{M}}_{i}[{\tilde{\mathbf{y}}}_{i}-\mathbf{P}(t){\alpha}_{i}(t)]\phantom{\rule{0ex}{0ex}}={[{\mathbf{P}}^{T}(t)\mathbf{W}(t){\mathbf{M}}_{i}\mathbf{P}(t)]}^{-1}{\mathbf{P}}^{T}(t)\mathbf{W}(t){\mathbf{M}}_{i}{\tilde{\mathbf{y}}}_{i},$$

(9)

where the polynomial matrix $\mathbf{P}(t)$ is given by $\mathbf{P}(t)={[{({t}_{i}-t)}^{j-1}]}_{i,j}$ for any $t$. The kernel weight matrix $\mathbf{W}$ is given by $\mathbf{W}(t)=\mathrm{diag}{\{{(Ns)}^{-1}K[({t}_{i}-t)/s]\}}_{i,i}$ where $K$ is a symmetric kernel function and $s$ is a smoothing parameter. We assume that the data are missing at random (MAR); i.e., $P[{m}_{i,t}=1|{y}_{i}(t)]=\mathfrak{p}(t)\in [\mathrm{0,1}]$.^{51} Then, the estimator is unbiased as $N\to \infty $ but has lower efficiency (in other words, the increase in $\mathrm{var}[{\widehat{f}}_{i}(t)|\mathbf{W},{\mathbf{M}}_{i}]$) as the observation probability $\mathfrak{p}(t)$ decreases.^{47}

1. Thomason M. E., et al. , “Cross-hemispheric functional connectivity in the human fetal brain,” Sci. Transl. Med.
5, 173ra24 (2013).http://dx.doi.org/10.1126/scitranslmed.3004978 [PMC free article] [PubMed]

2. Jakab A., et al. , “Fetal functional imaging portrays heterogeneous development of emerging human brain networks,” Front. Hum. Neurosci.
8, 1–17 (2014).http://dx.doi.org/10.3389/fnhum.2014.00852 [PMC free article] [PubMed]

3. Andescavage N. N., du Plessis A., Limperopoulos C., “Advanced MR imaging of the placenta: exploring the in utero placentabrain connection,” Semin. Perinatol.
39(2), 113–123 (2015).http://dx.doi.org/10.1053/j.semperi.2015.01.004 [PMC free article] [PubMed]

4. Buxton R. B., Introduction to Functional Magnetic Resonance Imaging: Principles and Techniques, 2nd ed., Cambridge University Press, Cambridge, UK: (2009).

5. Sørensen A., et al. , “Changes in human fetal oxygenation during maternal hyperoxia as estimated by BOLD MRI,” Prenat. Diagn.
33, 141–145 (2013).http://dx.doi.org/10.1002/pd.v33.2 [PubMed]

6. Sørensen A., et al. , “Changes in human placental oxygenation during maternal hyperoxia estimated by blood oxygen level-dependent magnetic resonance imaging (BOLD MRI),” Ultrasound Obstet. Gynecol.
42, 310–314 (2013).http://dx.doi.org/10.1002/uog.12395 [PubMed]

7. Malamateniou C., et al. , “Motion-compensation techniques in neonatal and fetal MR imaging,” Am. J. Neuroradiol.
34, 1124–1136 (2013).http://dx.doi.org/10.3174/ajnr.A3128 [PubMed]

8. Friston K. J., et al. , “Movement-related effects in fMRI time-series,” Magn. Reson. Med.
35, 346–355 (1996).http://dx.doi.org/10.1002/mrm.v35.3 [PubMed]

9. Hajnal J. V., et al. , “Artifacts due to stimulus correlated motion in functional imaging of the brain,” Magn. Reson. Med.
31, 283–291 (1994).http://dx.doi.org/10.1002/(ISSN)1522-2594 [PubMed]

10. Patel A. X., et al. , “A wavelet method for modeling and despiking motion artifacts from resting-state fMRI time series,” NeuroImage
95, 287–304 (2014).http://dx.doi.org/10.1016/j.neuroimage.2014.03.012 [PMC free article] [PubMed]

11. Yancey S. E., et al. , “Spin-history artifact during functional MRI: potential for adaptive correction,” Med. Phys.
38(8), 4634 (2011).http://dx.doi.org/10.1118/1.3583814 [PubMed]

12. Seshamani S., et al. , “A method for handling intensity inhomogenieties in fMRI sequences of moving anatomy of the early developing brain,” Med. Image Anal.
18, 285–300 (2013).http://dx.doi.org/10.1016/j.media.2013.10.011 [PMC free article] [PubMed]

13. Ferrazzi G., et al. , “Resting state fMRI in the moving fetus: a robust framework for motion, bias field and spin history correction,” Neuroimage
101, 555–568 (2014).http://dx.doi.org/10.1016/j.neuroimage.2014.06.074 [PubMed]

14. Sontag L., Wallace R., “The movement response of the human fetus to sound stimuli,” Child Development
6(4), 253–258 (1935).http://dx.doi.org/10.2307/1125396

15. Bekedam D., et al. , “The effects of maternal hyperoxia on fetal breathing movements, body movements and heart rate variation in growth retarded fetuses,” Early Human Development
27, 223–232 (1991).http://dx.doi.org/10.1016/0378-3782(91)90196-A [PubMed]

16. Cox R. W., “AFNI: software for analysis and visualization of functional magnetic resonance neuroimages,” Comput. Biomed. Res.
29, 162–173 (1996).http://dx.doi.org/10.1006/cbmr.1996.0014 [PubMed]

17. Woods R. P., et al. , “Automated image registration: I. General methods and intrasubject, intramodality validation,” J. Comput. Assisted Tomogr.
22, 139–152 (1998).http://dx.doi.org/10.1097/00004728-199801000-00027 [PubMed]

18. Woods R. P., et al. , “Automated image registration: II. Intersubject validation of linear and nonlinear models,” J. Comput. Assisted Tomogr.
22, 153–165 (1998).http://dx.doi.org/10.1097/00004728-199801000-00028 [PubMed]

19. Smith S., et al. , “FSL: new tools for functional and structural brain image analysis,” NeuroImage
13, 249 (2001).http://dx.doi.org/10.1016/S1053-8119(01)91592-7

20. Friston K. J., Jezzard P., Turner R., “Analysis of functional MRI time-series,” Hum. Brain Mapp.
1, 153–171 (1994).http://dx.doi.org/10.1002/hbm.v1:2

21. Friston K. J., et al. , “Spatial registration and normalization of images,” Hum. Brain Mapp.
3, 165–189 (1995).http://dx.doi.org/10.1002/hbm.v3:3

22. Glenn O. A., Barkovich A. J., “Magnetic resonance imaging of the fetal brain and spine: an increasingly important tool in prenatal diagnosis, part 1,” AJNR. Am. J. Neuroradiol.
27, 1604–1611 (2006). [PubMed]

23. D’Ercole C., et al. , “Prenatal diagnosis of fetal cerebral abnormalities by ultrasonography and magnetic resonance imaging,” Eur. J. Obstet. Gynecol. Reprod. Biol.
50, 177–184 (1993).http://dx.doi.org/10.1016/0028-2243(93)90198-L [PubMed]

24. Powell M. C., et al. , “Magnetic resonance imaging (MRI) in obstetrics. II. Fetal anatomy,” BJOG Int. J. Obstet. Gynaecol.
95, 38–46 (1988).http://dx.doi.org/10.1111/bjo.1988.95.issue-1 [PubMed]

25. Weinreb J. C., et al. , “Magnetic resonance imaging in obstetric diagnosis,” Radiology
154, 157–161 (1985).http://dx.doi.org/10.1148/radiology.154.1.3880601 [PubMed]

26. Gowland P., Fulford J., “Initial experiences of performing fetal fMRI,” Exp. Neurol.
190(Suppl. 1), 22–27 (2004).http://dx.doi.org/10.1016/j.expneurol.2004.06.022 [PubMed]

27. Schöpf V., et al. , “Watching the fetal brain at ‘rest’,” Int. J. Dev. Neurosci.
30(1), 11–17 (2012).http://dx.doi.org/10.1016/j.ijdevneu.2011.10.006 [PubMed]

28. Fulford J., et al. , “Fetal brain activity in response to a visual stimulus,” Hum. Brain Mapp.
20(4), 239–245 (2003).http://dx.doi.org/10.1002/(ISSN)1097-0193 [PubMed]

29. Fulford J., et al. , “Fetal brain activity and hemodynamic response to a vibroacoustic stimulus,” Hum. Brain Mapp.
22(2), 116–121 (2004).http://dx.doi.org/10.1002/(ISSN)1097-0193 [PubMed]

30. Thomason M. E., et al. , “Age-related increases in long-range connectivity in fetal functional neural connectivity networks in utero,” Dev. Cogn. Neurosci.
1–9 (2014).http://dx.doi.org/10.1016/j.dcn.2014.09.001 [PMC free article] [PubMed]

31. You W., et al. , “Robust motion correction and outlier rejection of in vivo functional MR images of the fetal brain and placenta during maternal hyperoxia,” Proc. SPIE
9417, 941700 (2015).http://dx.doi.org/10.1117/12.2082451 [PMC free article] [PubMed]

32. Juntu J., et al. , “Bias field correction for MRI images,” in Computer Recognition Systems, Kurzyski M., et al., editors. , Eds., Vol. 30, pp. 543– 551, Advances in Soft Computing, Chapter Part IV, Springer Berlin Heidelberg, Berlin, Heidelberg: (2005).

33. Tustison N. J., et al. , “N4ITK: improved N3 bias correction,” IEEE Trans. Med. Imaging
29, 1310–1320 (2010).http://dx.doi.org/10.1109/TMI.2010.2046908 [PMC free article] [PubMed]

34. Jenkinson M., et al. , “Improved optimization for the robust and accurate linear registration and motion correction of brain images,” NeuroImage
17, 825–841 (2002).http://dx.doi.org/10.1006/nimg.2002.1132 [PubMed]

35. Yushkevich P. A., et al. , “User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability,” NeuroImage
31(3), 1116–1128 (2006).http://dx.doi.org/10.1016/j.neuroimage.2006.01.015 [PubMed]

36. Crum W. R., Hartkens T., Hill D. L. G., “Non-rigid image registration: theory and practice,” Br. J. Radiol.
77, S140–S153 (2004).http://dx.doi.org/10.1259/bjr/25329214 [PubMed]

37. Szeliski R., “Image alignment and stitching: a tutorial,” Found. Trends. Comp. Graphics Vision
2(1), 1–104 (2006).http://dx.doi.org/10.1561/0600000009

38. Serra J. P., Cressie N. A. C., “Image analysis and mathematical morphology,” in Image Analysis and Mathematical Morphology, Vol. 1, Academic Press, New York, NY: (1982).

39. Summers P. E., et al. , “A quantitative comparison of BOLD fMRI responses to noxious and innocuous stimuli in the human spinal cord,” NeuroImage
50, 1408–1415 (2010).http://dx.doi.org/10.1016/j.neuroimage.2010.01.043 [PubMed]

40. Keehn B., et al. , “Functional connectivity for an “island of sparing” in autism spectrum disorder: an fMRI study of visual search,” Hum. Brain Mapp.
34, 2524–2537 (2013).http://dx.doi.org/10.1002/hbm.v34.10 [PMC free article] [PubMed]

41. Becerra L., et al. , “Robust reproducible resting state networks in the awake rodent brain,” PLoS One
6, e25701 (2011).http://dx.doi.org/10.1371/journal.pone.0025701 [PMC free article] [PubMed]

42. Ahn H. K., Brass P., Shin C. S., “Maximum overlap and minimum convex hull of two convex polyhedra under translations,” Comput. Geom. Theory Appl.
40, 171–177 (2008).http://dx.doi.org/10.1016/j.comgeo.2007.08.001

43. Ahn H.-K., et al. , “Overlap of convex polytopes under rigid motion,” Comput. Geom.
47, 15–24 (2014).http://dx.doi.org/10.1016/j.comgeo.2013.08.001

44. Grimmett G., Stirzaker D., Probability and Random Processes, 3rd ed., Oxford University Press, Oxford, UK: (2001).

45. Papoulis A., Pillai S., Probability, Random Variables and Stochastic Processes, 3rd ed., McGraw-Hill Education (India) Pvt Ltd., New York, NY: (2002).

46. Rosner B., “Percentage points for a generalized ESD many-outlier procedure,” Technometrics
25, 165–172 (1983).http://dx.doi.org/10.1080/00401706.1983.10487848

47. Pérez-González A., Vilar-Fernández J. M., González-Manteiga W., “Asymptotic properties of local polynomial regression with missing data and correlated errors,” Ann. Inst. Stat. Math.
61, 85–109 (2007).http://dx.doi.org/10.1007/s10463-007-0136-2

48. Abaci Turk E., et al. , “Automated ROI extraction of placental and fetal regions for 30 minutes of EPI BOLD acquisition with different maternal oxygenation episodes,” Proc. Intl. Soc. Mag. Reson. Med.
23, 0639 (2015).

49. Andersson J. L. R., Jenkinson M., Smith S., “Non-linear registration aka spatial normalisation,” Technical Report, FMRIB Centre, Oxford, United Kingdom: (2007).

50. Field A. S., et al. , “False cerebral activation on BOLD functional MR images: study of low-amplitude motion weakly correlated to stimulus,” Am. J. Neuroradiol.
21(8), 1388–1396 (2000). [PubMed]

51. Acock A. C., “Working with missing values,” J. Marriage Fam.
67, 1012–1028 (2005).http://dx.doi.org/10.1111/jomf.2005.67.issue-4

Articles from Journal of Medical Imaging are provided here courtesy of **Society of Photo-Optical Instrumentation Engineers**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |