Blood oxygenation level-dependent (BOLD) contrast is a $T^{*}_{2}$-weighted sequence. It is affected by variations in concentration of vascular oxygentation in the blood volume and magentic field inhomogeneities. Quantifying $T_{2}^{*}$ enables the determination of oxygen saturation by leveraging the relationationship between $T_{2}^{*}$ and deoxyhemoglobin (Sinding et al. (2016, 2017)). In FGR pregnancies, the placenta is hypoxic, displaying a reduced $T_{2}^{*}$ value which can be used as an FGR biomarker (Robinson et al. (1998); Jiang et al. (2013)). Despite the potential of BOLD-MRI in measuring oxygen saturation, its use has not yet been validated in diagnosis of FGR and interpretation of the placental BOLD signal is complicated by several factors that influence changes in the signal (Sinding et al. (2018); Uğurbil et al. (2000); Chalouhi and Salomon (2014); Sørensen et al. (2015); Turk et al. (2020)). Considering this, and due to the requirements of a gradient echo acquisition, $T_{2}^{*}$ relaxometry is not quantified in the current research.

Instead, T${}_{2}$ relaxometry provides structural, functional, and morphological tissue information as T${}_{2}$ transverse relaxation times depend on several factors encompassing water binding, macromolecular concentration, and most importantly, blood oxygenation levels (Derwig et al. (2013); Saini et al. (2020)). Previous literature has shown that the placental T${}_{2}$ times in SGA or FGR pregnancies are reduced with respect to normal pregnancies (Derwig et al. (2013)). T${}_{2}$ relaxation times have been used to assess placental function in various applications (Melbourne et al. (2019, 2016a); Jacquier and Salomon (2021); Stout et al. (2021)).

This study extended on previous placental research by producing T${}_{2}$ maps for different fetal organs. It was hypothesised that because of the brain-sparing effect, certain organs would have lower oxygen levels in FGR compared to healthy pregnancies, and thus reduced T${}_{2}$ measurements would be extracted from blood flowing through non-prioritised organs in FGR pregnancies. Portnoy et al. demonstrated the precise relationship between blood T${}_{2}$ relaxation times and oxygen saturation (Portnoy et al. (2017)) by making use of the Luz-Meiboom model, given by Equation 1, presenting the exponential relationship,

where $R_{2}=\frac{1}{T_{2}}$, $R_{2,ery}$ is the erythrocyte (red blood cell) relaxation rate that depends on oxygen saturation, $Hct$ is the hematocrit (proportion of red blood cells in blood), and $R_{2,plas}$ is the plasma relaxation rate.

Diffusion-weighted (DW) MRI is a valuable method for investigating the fetal brain-sparing effect; providing measures of brain maturation and detection of brain lesions (Arthurs et al. (2017)). This is attained by measuring water diffusion, which yields corresponding apparent diffusion coefficient (ADC) values. Arthurs et al. established differences between healthy and severe FGR fetuses, frequently leading to the clinical decision of early delivery induction in the latter group (Arthurs et al. (2017)). The time between the MRI examination and delivery was, on average, 7.69 weeks earlier for the FGR group compared to the healthy, thus highlighting the potential of DW-MRI for accurate diagnosis of growth restricted cases - allowing for appropriate management plans to be put in place.

Dynamic contrast-enhanced (DCE) MRI can spatially and quantitatively characterise maternal perfusion of placental insufficiency and tissue vasculature (Ingram et al. (2018); Schrauben et al. (2019); Frias et al. (2015)). It describes the delivery of contrast agent to the maternal side and its transfer into the fetal blood pool in order to distinguish between individual vascular units of the placenta. DCE-MRI is the current gold standard for quantitative descriptions of vascular function (Frias et al. (2015); Schabel et al. (2016)). Nonetheless, this technique has significant drawbacks as it requires an exogenous contrast. The clearance of contrast from the feto-placental system still requires further research. To that end, an imaging technique which does not include any safety concerns for the mother and fetus is more pertinent.

Multi-compartment models refer to advanced mathematical models that separate the signal contributions from different tissue types (Aughwane et al. (2020b)). Diffusion-relaxation models are growing in popularity and have found multiple applications such as in neuroimaging (Kim et al. (2017)), and more recently in placental imaging, encompassing the assessment of placental function in FGR (Melbourne et al. (2019, 2016a); Hutter et al. (2019); Jacquier and Salomon (2021)). A thorough overview of these techniques is provided in (Slator et al. (2021)).

The DECIDE model identifies and separates the $T_{2}$ values corresponding to the fetal and maternal blood, enabling the quantification of fetal blood oxygen saturation. The precise mechanisms and assumptions describing the DECIDE and the Extended Intravoxel Incoherent Motion Model (IVIM) models are discussed in Section 3.2.

Supervised Machine Learning (ML) refers to the employment of a predictive model with an assumed relationship between the input (features) and output (labels) variables. Its prominence in medical imaging has been significantly established in recent years, particularly in computer-aided diagnosis (Erickson et al. (2017); Giger (2018)), due to the rise of ‘Big Data’ and available computer power. Its contribution to “intelligent imaging” is by virtue of its potential to advance and enhance detection and diagnosis of complex disorders, risk assessment, and therapy response (Schoepf et al. (2007); Dundar et al. (2008); Summers (2010); Mitchell et al. (2008)). The advantages of ML stem from its ability to draw connections and identify patterns between variables, surpassing human perception. However, its attribute as a ‘black-box function’ makes it difficult to interpret the results from ML models and determine how features are used to arrive at predictions, thus ensuing in a lack of clinician trustworthiness in the models. Nonetheless, ML can be leveraged to assimilate information from datasets where the relationship between the input and output variables are unknown and to select the best features for a certain prediction. It can be used for decision support by aiding clinicians in interpreting medical imaging findings rather than relying entirely on the model predictions alone.

ML enables the consolidation and unravelling of complex biomedical and healthcare data that overcomes the limitations of traditional statistical methods. The aim of these algorithms is to provide solutions to clinical problems by learning statistical associations of the features extracted from the images (Shen et al. (2017)).

Current screening and diagnostic tools for FGR remain suboptimal (Audette and Kingdom (2018)). Delivery of improved clinical outcomes requires greater understanding of the multifactorial pathogenesis in early-onset FGR and distinguishing features or biomarkers of the condition (Audette and Kingdom (2018)). Analysis of a combination of multiple FGR indicators (Gordijn et al. (2018, 2016b); Beune et al. (2018)), can be achieved through use of ML methods. Supervised ML models are increasingly being employed for early prediction and diagnosis of pregnancy conditions, including intrauterine growth restriction, pre-eclampsia, risk of stillbirth, preterm pregnancy, and gestational diabetes (Crockart et al. (2021); Burgos-Artizzu et al. (2020); Caly et al. (2021); Khatibi et al. (2021); Marić et al. (2020); Ye et al. (2020); Koivu and Sairanen (2020)).

Recent work conducted by (Arabi Belaghi et al. (2021)) compared the performance of logistic regression and artificial neural networks in predicting overall and spontaneous preterm birth, on a dataset of 112,963 nulliparous women (singleton gestation) who delivered between 20-42 weeks gestation. The predictors included socio-demographic variables correlated with the risk of preterm birth, such as maternal age, income, education, race, folic acid use, etc. The prediction accuracy of both models in the first trimester was ambiguous. But by incorporating complications during pregnancy as additional predictors, the authors established a 20% increase in the area under the curve (AUC) from the receiver operating characteristic curve (ROC) for artificial neural networks in the validation sample compared to the logistic regressor during the second trimester (80% vs. 60%). The prediction performance of this work cannot be directly compared to our study, given the substantial difference in sample size (being several orders of magnitude smaller), which greatly influences the statistical power of the study. Therefore, our study should be viewed only as a preliminary study, as gaining concrete and detailed information regarding model performance and generalisability, and the features driving each ML model cannot be easily extracted as in (Arabi Belaghi et al. (2021)), where the statistical methods were applied to a much larger dataset.

Research into the prediction of stillbirth by (Yerlikaya et al. (2016)) employed a multivariate logistic regression analysis to deduce the contributions of varying maternal characteristics and medical history in stillbirth prediction. Correspondingly, Trudell et al. generated models for the prediction of stilbirth using backward stepwise logistic regression (Trudell et al. (2017)). Both groups leveraged highly similar maternal demographics and medical history and concluded similar prediction performances ranging between 64% to 67% AUC.

Despite the comparable performance of conventional logistic regression and ML methods for diagnosis predictions in previous literature (Yerlikaya et al. (2016); Trudell et al. (2017); Koivu and Sairanen (2020); Ye et al. (2020)), the former assumes linearity and independence between the features. As such, we extended our previous methods (Zeidan et al. (2021)) which implemented logistic regression to diagnose FGR and assess its severity, to a convolutional neural network (CNN). Deep learning algorithms draw on higher-level features extracted from the lower-level features of input data (Bengio (2012)). In particular, the benefits can be observed in supervised learning due to the scalability of deep neural networks and feature learning abilities. The use of a CNN allows us to explore both spatial and intensity relationships at a voxel-wise level for each of our parameter maps - information which is otherwise excluded when employing simple logistic regression models over averaged maps. Thus, we aim to maximise feature extraction from our parameter maps via a CNN, where feature representation is more accurate to the underlying maps for each organ, as no high-level averaging takes place. Nonetheless, it is important to acknowledge that a considerable amount of data is crucial to obtaining robust ML models.

Patient MRI scans of voxel resolution 1.9x1.9x6$mm$ were acquired using the acquisition parameters from (Melbourne et al. (2019)), where b-values and echo times are varied in pairs (enabling both T${}_{2}$ relaxometry and DW-MRI fitting), using a 1.5 T Siemens Avanto and performed under free-breathing. The dataset consisted of 12 early-onset FGR (Gordijn et al. (2016a)) ranging between [24${}^{+2}$, 33${}^{+6}$] gestation weeks${}^{+days}$, and 12 control pregnancies with MR data ranged between [25${}^{+1}$, 34${}^{+0}$] GA interval, (Median 28${}^{+4}$wks$\pm$3${}^{+3}$wks) respectively. Specific details on subject inclusion criteria are available in (Aughwane et al. (2020b)). The study was approved by the UK National Research Ethics Service and all participants gave written informed consent (REC reference 15/LO/1488).

There are biological mechanisms that may cause differences in the distribution of blood perfusion throughout the fetus in FGR. To investigate this, manual segmentation of the placenta, liver, lungs and brain was accomplished using the open-source ITK-SNAP application (image segmentation). The resultant 3D mask files were used within the NiftyFit package (Melbourne et al. (2016b)) for multi-parametric model-fitting (Melbourne et al. (2019)), and to perform texture analysis.

Model fitting techniques were applied to each organ segmentation over the averaged region of interest (ROI) signal and on a voxelwise scale, yielding quantitative metrics for both approaches. Non-linear least squares were used to perform the fitting, with voxelwise fitting being initialised with the ROI parameter estimates - enhancing signal-to-noise ratio (SNR) by reducing the changes of fitting to local minima. A range of models were explored, including simple T${}_{2}$ and ADC estimation, as well as more complex models based on IVIM (Le Bihan et al. (1986)) and DECIDE (Melbourne et al. (2019)). Investigated in this research were parameters linked to diffusion, but they do not represent diffusion directly. The simplest T${}_{2}$ model fitting describes the MRI signal as

where $TE$ are the echo times, $S$ is the measured signal, and $S_{0}$ the baseline signal. Regarding simple ADC fitting, this is accomplished using

where b are the b-values. Thus, the acquired data requires varying $TE$ and b-values to allow for dual ADC and T2 model fitting.

The IVIM model describes perfusion as a pseudodiffusion process (represented by a pseudodiffusion coefficient, D${}^{*}$), by characterising the collective motion of blood water molecules within the vessel network as a random walk. The IVIM model also incorporates “true” diffusion of water molecules (ADC), modelling the signal as

where $f$ is the perfusion fraction (volume occupied by incoherently flowing blood in a given voxel) and $b$ is the b-value (Le Bihan (2019)). We refer to this model as Standard IVIM (Eq. 4). This can be extended to incorporate T${}_{2}$ relaxometry as

We refer to this model (Eq. 5) as T2 IVIM. However, this model presents inherent limitations, as it assumes both vascular and tissue compartments (parametrised by pseudo-diffusion and true diffusion coefficients) have the same T${}_{2}$ value, leading to an overestimation of the pseudo-diffusion volume fraction $f$ with increasing echo time ($t$) (Jerome et al. (2016)). Thus, the presented analysis incorporates more complex models, accounting for varying blood and tissue T${}_{2}$ values:

with $f$ being the perfusion fraction, T${}_{2p}$ and T${}_{2t}$ being the transverse relaxation time for the pseudo-diffusion compartment (blood) and true diffusion compartment (tissue), respectively (Jerome et al. (2016)). We refer to this model as Extended $2\times$T2 IVIM (Eq. 6).

The DECIDE model (Melbourne et al. (2019)) was also applied specifically to the placenta, which assumes three compartments with distinct diffusivity and relaxivity: fetal capillaries, trophoblast space and maternal blood pool. This model, given by Equation 7, enables computation of novel placental biomarkers including maternal fetal blood volume ratio and fetal blood saturation.

Here, T${}_{2}^{fb}$, T${}_{2}^{mb}$ and T${}_{2}^{ts}$ represent the transverse relaxation times for fetal blood, maternal blood and trophoblast space, respectively; and $\nu$ is the maternal blood volume fraction. R${}_{2}^{mb}$ and R${}_{2}^{ts}$ are fixed known, $(240ms)^{-1}$ and $(46ms)^{-1}$ respectively at 1.5T), taken from (Melbourne et al. (2019)).

The aim of texture analysis was to examine the spatial arrangement of intensities in the segmented organs using in-house software developed in MATLAB (The MathWorks Inc., Natick, MA). To perform the texture analysis, a grey level co-occurrence matrix (GLCM) was computed to provide insight into the spatial interaction of neighbouring pixels. Haralick features are statistical features extracted from the GLCM to describe the overall image texture using measures encompassing energy, entropy, correlation, contrast, variance, and homogeneity (Haralick et al. (1973)):

Energy: This measure is extracted from the angular second moment, which calculates the grey level local uniformity,

where $i$ and $j$ represent the image dimensions, and $p_{d}(i,j)$ corresponds to an element of the normalised GLCM.

Entropy: A statistical measure of randomness.

Correlation: A measurement of the similarity between neighbouring pixels,

where $\mu_{x}$; $\mu_{y}$ are the means and $\sigma_{x}$; $\sigma_{y}$ are the standard deviations.

Contrast: The number of grey levels that exist in the scan.

Variance: A measure of variability.

Homogeneity: The number of changes of intensity that appear in a region of interest.

We hypothesised that these six Haralick features could be used to discern between FGR and appropriately grown fetuses due to a lower SNR present in FGR fetuses as a result of lower $T_{2}$ and decreased oxygen saturation (Portnoy et al. (2017)). We expected that the lower signal intensities in FGR compared to controls would be especially evident in the placenta and fetal liver and correlate directly with placental insufficiency (Aughwane et al. (2020a); Kessler et al. (2009)). For instance, this would be reflected in the computed Haralick features by observing lower contrast values in FGR fetuses in comparison to the controls. Decreased contrast in the ROI would equate to an increase in homogeneity.

These features were computed for each subject on the most significant parameter maps for each organ (as determined by the t-tests described in Section 3.4 with a p-value cut-off of 0.05), as well as the b=0 volume with lowest echo time from the original IVIM T${}_{2}$-weighted MRI scan; this yielded interpretable texture descriptors (Haralick et al. (1973); Bharati et al. (2004)). The images were quantised into grey level bins of fixed equal width for between-subject texture feature value comparisons. Single-factor analysis of each feature was conducted between the FGR and control patients. Results from the texture analysis were then concatenated by considering the mean and max of each Haralick feature.

The model fitting maps provide voxelwise information for each of the parameters optimised for. We simplified this information by considering the mean, max, variance and mode of each of voxelwise map. This yielded reduced parameters to be used for subsequent classical ML-based assessments.

We performed statistical analysis on these simplified model fitting parameters and on the Haralick features, in order to identify the most significant features in differentiating between the control and FGR cohorts. A Shapiro–Wilk test was used to confirm normality of the parameters obtained from the model fitting on a patient-by-patient basis to verify it was justifiable to run a t-test on them. In particular, the Shapiro-Wilk test was selected for its efficacy on small sample sizes. The test was run on the distribution of the model fitted parameters for each of the organ ROIs done over all of the samples split between the two cohorts.

T-tests were then carried out between the two cohorts for all the model fitted parameters, Haralick features, and organ ratio parameters. Results with p-value less than 0.05 indicated statistically significant differences between the control and FGR group means.

We used these significant parameters for training simple classical machine learning models on classification (control or FGR) and regression (GA at birth, time from scan until birth and baby weight), as detailed in the following sections.

Following these statistical tests, we aimed to explore the use of these significant features (p-value $<0.05$ in distinguishing between controls and FGR) as potential FGR biomarkers for severity assessments. For this, we conducted various ML training experiments, employing a binary classifier for diagnosis prediction (control or FGR), and simple regressors to predict GA at birth, time from scan until birth, and baby weight.

Our training experiments explored the most appropriate use of our data to achieve optimal results. For this, we trained each model first using exclusively model fitting data (mean, max, variance and mode of each of voxelwise map), followed by exclusive training using Haralick features, and finally combining both model fitting data and Haralick features. Only the features with a p-value$<0.05$ in differentiating between controls and FGR cohorts were employed.

We employed logistic regression for binary classification, using a stochastic average gradient (SGA) solver that supports the L1 regularisation to minimise the cross-entropy loss function.

This classifier models the conditional probability of an FGR or non-FGR (control) diagnosis, Y, given input features, X (model fitting data and Haralick features), by applying a sigmoid function to the output of a decision function $h(\textbf{x})=\textbf{w}^{T}\textbf{X}$, which ensures an output between 0 and 1:

where X is the input feature vector, and w is the learnt weight vector. These probability scores (i.e. the output for Eq. 14) are mapped to discrete classes with a decision boundary of 0.5, that is, an output probability$<0.5$ indicates an FGR diagnosis, while an output probability$\geq 0.5$ specifies a non-FGR diagnosis.

The optimal regularisation parameters (found via a grid search) were an L1 ratio of 0, i.e. L2 regularisation for all classifiers; and a regularisation strength ($C$) of 0.001 for the classifier trained exclusively on model fitting features, as well as the joint model (Haralick and model fitting features), while the model trained only on Haralick features yielded a $C=0.25$.

Based on RFECV, we used 44 out of 84 features for the classifier trained on model fitting data; 34 out of 53 features for the classifier trained on Haralick features; and 118 out of 137 for the classifier trained on both feature types.

We trained three multi-variate linear regressors to predict GA at delivery, time interval between scan and delivery, and baby weight, as these variables ($\hat{\textbf{y}}$) are potential indicators of FGR severity. Thus we fitted a linear equation $\hat{\textbf{y}}=\textbf{X}\textbf{w}$ to our feature matrix X, minimising the sum of squared errors between predicted and expected target values (including L1 and L2 regularisation) in order to find the weights, w.

Refer to Table 3 in Section 4.5 for information regarding model tuned hyperparameters and number of selected features for each regressor.

The regression methods described in Section 3.7, use data which statistically shows differences between FGR and healthy (p-value$<$0.05), followed by RFECV feature selection, to further reduce the noise present. However, the features used (in Section 3.7), particularly the model fitting features, drastically reduce the amount of parameter maps information: by taking single statistical values over whole voxelwise maps (i.e. mean, max, min, mode), important spatial relationships and detailed voxel-level information may be eliminated. The Haralick features do contain information regarding spatial arrangements and intensity relationships, which supports our previous method.

In an attempt to make use of this detailed information contained within each parameter map, we explored the potential of a Convolutional Neural Network (CNN) for severity assessment, aiming to predict the same regression variables as with our simple ML models (GA at birth, time interval from scan until delivery, and baby weight).

Figure 2 depicts examples of the parameter maps obtained from the model fitting techniques. The lower parameter map intensities in FGR compared to that in the controls is indicative of hypoperfusion and low oxygen saturation levels in these fetal organs. The T${}_{2}$ maps display pronounced differences in the signal intensities of both cohorts.

The most significant parameters in identifying differences between controls and FGR fetuses were the perfusion fraction, S${}_{0}$, pseudo-diffusion coefficient (D${}^{*}$), and T${}_{2}$ as given in Table 1. The placenta and liver were determined to be the most influential organs in diagnosing FGR.

The results for the parameter feature importances in Table 1, specify that there were no significant differences detectable in the fetal brain and lungs between normal and FGR fetuses, especially compared to the placenta and liver, where differences were significant.

Evaluation of the resulting Haralick features corroborated the degree of effect on the placenta in FGR, particularly using the Extended T${}_{2}$ IVIM map and its mean variance. The brain was the least significantly different organ in this analysis. Greater mean variance in the signal from the Extended T${}_{2}$ IVIM model of the healthy cohort (refer to Figure 4(a)), is indicative of increased heterogeneity in FGR placentas. The max correlation of the liver perfusion fraction in the controls in Figure 4(d) reflects larger intensity differences compared to FGR. This is a significant feature to consider in the Standard IVIM model when studying the liver in FGR, especially given that the notches do not overlap between the cohorts.

Referring to the results presented in Table 2, the classifier performs best when trained exclusively on model fitting data, achieving a prediction accuracy of 100% in testing, and thus a precision and recall score of 1.0.

This is further validated by the cross validated accuracy on training set, with a standard deviation of only 10% across folds, hinting at optimal model generalisability.

Given our optimal test set classification results (see Table 2), we qualitatively assess the most important features driving each classifier model. These were obtained via Recursive Feature Elimination (RFE). We obtained the exact same top five features for both the logistic regressor trained exclusively on model fitting data, and the logistic regressor trained on both Haralick features and model fitting data.

Figure 5b shows distinct differences between controls and FGR cohorts. Here, control subjects display a much higher placenta and liver perfusion relative to the lungs, compared to FGR subjects. This is an indicator that in FGR, both the placenta and the lungs are much less perfused than other vital organs such as the lungs.

The differences between the performance of the model fitting techniques can be inferred from Figure 5c. It showcases the ‘liver/lungs median f’ plotting the Extended 2xT2 IVIM model against the Standard IVIM model. When projected onto each axis, the x-axis (i.e. the Standard IVIM model), permits a more accurate linear classification between the cohort compared to the Extended 2xT2 IVIM model. Seemingly, the Standard IVIM model produces a better fit to our data (for this particular parameter) than the more complex models with additional parameters, potentially due to the noise present.

Table 3 includes our test set and cross validated regressor results. In accordance with our classifier results, the models with highest performance are those trained on model fitting features, excepting predictions for time from scan until delivery, where the combined model displays an insignificantly lower root mean square error (RMSE) on test set compared to the model trained exclusively on model fitting data.

Figure 6 depicts test set regression predictions for our best performing regressors, against true labels. Qualitatively, the test set predictions that mostly resemble the true data points are from the time interval between MRI and delivery (Figure 6b), however this has two important outliers. It is complex to comment on the significance of this, given the extremely small test set. The plot depicting baby weight predictions (Figure 6c) visually appears as the worst fit, however the value range for this variable is much larger, which may partially explain this. Additionally, the most important outliers for baby weight are predictions which are lower than the actual baby weight, which is clinically significant: it is best to overestimate the severity than underestimate it.

The ResNet results for severity assessment are included in Table 4. Both GA at birth and time interval between scan and delivery resulted in a much higher RMSE than those obtained from our simple classical linear regressors (see Section 4.5).