# RECENT ADVANCES IN DOPPLER SIGNAL PROCESSING AND MODELLING TECHNIQUES FOR FETAL MONITORING

EDITED BY : Ahsan H. Khandoker, Faezeh Marzbanrad and Yoshitaka Kimura PUBLISHED IN : Frontiers in Physiology

#### Frontiers Copyright Statement

© Copyright 2007-2018 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88945-536-2 DOI 10.3389/978-2-88945-536-2

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# RECENT ADVANCES IN DOPPLER SIGNAL PROCESSING AND MODELLING TECHNIQUES FOR FETAL MONITORING

Topic Editors:

Ahsan H. Khandoker, Khalifa University, United Arab Emirates Faezeh Marzbanrad, Monash University, Australia Yoshitaka Kimura, Tohoku University, Japan

Image: Dmitry Naumov/Shutterstock.com

The guest editors of this eBook have accepted 10 very high-quality submissions for inclusion in a special issue of Frontiers in Physiology. The key difference between this eBook and contemporary fetal physiology related literature is that this Research Topic summarizes additional insights into the physiological link between physiologically understandable mathematical indices of fetal signals and the developing cardiovascular functions in fetal health and compromises. This book should be of considerable help to researchers, professionals in fetal monitoring device industries, academics, and graduate students from a wide range of disciplines. The text provides a comprehensive account of recent research in this emerging field and we anticipate that the concepts presented here will generate further research in this field.

Citation: Khandoker, A. H., Marzbanrad, F., Kimura, Y., eds. (2018). Recent Advances in Doppler Signal Processing and Modelling Techniques for Fetal Monitoring. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-536-2

# Table of Contents

*04 Editorial: Recent Advances in Doppler Signal Processing and Modeling Techniques for Fetal Monitoring* Ahsan H. Khandoker, Faezeh Marzbanrad and Yoshitaka Kimura

### 1. FETAL CARDIOGENIC SIGNALS


Saeed Abdulrahman Alnuaimi, Shihab Jimaa and Ahsan H. Khandoker

*17 Template-Based Quality Assessment of the Doppler Ultrasound Signal for Fetal Monitoring*

Camilo E. Valderrama, Faezeh Marzbanrad, Lisa Stroux and Gari D. Clifford

*27 Fetal Heart Sounds Detection Using Wavelet Transform and Fractal Dimension*

Elisavet Koutsiana, Leontios J. Hadjileontiadis, Ioanna Chouvarda and Ahsan H. Khandoker


Clarissa L. Velayo, Kiyoe Funamoto, Joyceline Noemi I. Silao, Yoshitaka Kimura and Kypros Nicolaides

# 2. FETAL HEART RATES AND VALVE TIMING INTERVALS

*57 A Hybrid EMD-Kurtosis Method for Estimating Fetal Heart Rate From Continuous Doppler Signals*

Haitham M. Al-Angari, Yoshitaka Kimura, Leontios J. Hadjileontiadis and Ahsan H. Khandoker


# 3. ANIMAL MODEL IN PERINATAL MONITORING

*92 Ultrasound Imaging of Mouse Fetal Intracranial Hemorrhage Due to Ischemia/Reperfusion*

Kenichi Funamoto, Takuya Ito, Kiyoe Funamoto, Clarissa L. Velayo and Yoshitaka Kimura

# Editorial: Recent Advances in Doppler Signal Processing and Modeling Techniques for Fetal Monitoring

#### Ahsan H. Khandoker <sup>1</sup> \*, Faezeh Marzbanrad<sup>2</sup> and Yoshitaka Kimura<sup>3</sup>

*<sup>1</sup> Department of Biomedical Engineering, Khalifa University, Abu Dhabi, United Arab Emirates, <sup>2</sup> Department of Electrical and Computer Engineering, Monash University, Clayton, VIC, Australia, <sup>3</sup> Graduate School of Medicine, Tohoku University, Sendai, Japan*

Keywords: fetal Doppler ultrasound, fetal electrocardiogram (FECG), fetal phonocardiogram, fetal monitoring, fetal screening

#### **Editorial on the Research Topic**

#### **Recent Advances in Doppler Signal Processing and Modeling Techniques for Fetal Monitoring**

The intersection of perinatal medicine and biomedical engineering is an emerging scientific research area. Understanding pathophysiological process of fetal development and improving the technologies used in perinatal diagnosis and intervention in clinical settings are vital for monitoring fetal well-being reliably. The major challenges in fetal Doppler, acoustic and electrical signal processing techniques which are often used in obstetrical instrumentations are poor specificity with high false positive rates and strong non-stationarities in abdominal derived signals. Therefore, more research is needed to explore the untapped potentials of abdominal sensor or lead based fetal signal (Doppler, ECG, Phonogram etc.) analyses and modeling.

In this research topic, selected emerging technologies in abdominal derived fetal signal processing for screening of fetal development and well-being are published. These papers include reviews, commentaries and technical contributions on: (1) challenges in fetal cardiac signal processing techniques; (2) nonlinear signal processing techniques in fetal Doppler signals; (3) novel techniques for detection of fetal heart sounds; (4) quality assessment of Doppler signals; (5) abdominal fetal electrocardiography technique and its application on growth restricted fetuses; (6) comparative study on fetal heart rates; (7) application of fetal Doppler ultrasound and fetal electrocardiography in pregnant animal model (mouse).

The guest editors of this research topic have accepted 10 very high-quality submissions for inclusion in this special issue. The key difference between this issue and contemporary fetal physiology related literature is that this research topic summarizes additional insights into the physiological link between physiologically understandable mathematical indices of fetal signals and the developing cardiovascular functions in fetal health and compromises. The summary of the research papers published in this topic is given below within three main categories.

### FETAL CARDIOGENIC SIGNALS

Cardiogenic such as mechanical, auscultation and electrical activities of fetal heart are acquired by abdominal Doppler, Phonogram and ECG technologies.

Frasch et al. (2017) in a general commentary shared the experience from a randomized controlled trial with computerized interpretation of fetal heart rate during labor (INFANT)

Edited and reviewed by: *Zbigniew R. Struzik, The University of Tokyo, Japan*

> \*Correspondence: *Ahsan H. Khandoker ahsank@ieee.org*

#### Specialty section:

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology*

Received: *22 March 2018* Accepted: *17 May 2018* Published: *05 June 2018*

#### Citation:

*Khandoker AH, Marzbanrad F and Kimura Y (2018) Editorial: Recent Advances in Doppler Signal Processing and Modeling Techniques for Fetal Monitoring. Front. Physiol. 9:691. doi: 10.3389/fphys.2018.00691*

**4**

collaborative group. They proposed the importance of new modalities of direct acquisition of predictive information about the fetal brain and heart so that a variety of fetal compromises could be prevented.

Alnuaimi et al. reviewed the most recent progress in fetal cardiac Doppler signal processing research. The review focuses on the shortcomings and advantages, which helps in understanding fetal Doppler cardiogram signal processing methods, and the related Doppler signal analysis procedures, by providing valuable clinical information, with a set of recommendations for future research directions.

One of the challenges of fetal monitoring techniques, such as Doppler ultrasound, is the susceptibility to noise affecting the signal quality, and subsequently the accuracy and reliability of the derived metrics (Stroux and Clifford, 2014; Marzbanrad et al., 2015; Valderrama et al., 2018). Valderrama et al. proposed an interesting study on how to assess the quality of Doppler signals, which are inherently non-stationary, and highly movement sensitive signals. The proposed template based method introduced a new signal quality index, which is validated on recordings from a low cost Doppler probe. This novel method might be useful in resource-poor settings or when operated by non-skilled operators.

Koutsiana et al. demonstrated a Wavelet Transform-based method combined with the Fractal Dimension to detect fetal heart sounds from abdominal phonogram signals. Fetal phonography is popular with low-cost option in many developing countries in the world.

Li et al. (2017) demonstrated an efficient method of extracting fetal and maternal ECG signals from two channels, which are attached on abdominal leads. The proposed method was then validated on a non-invasive fECG database in PhysioNet. The proposed technique might be useful in ambulatory monitoring of fetal vital signals at home.

Velayo et al. then showed a real application of abdominal lead based fetal ECG in growth-restricted fetuses. The results of the study confirmed that both QT and QTc intervals were significantly prolonged in growth-restricted fetuses as compared to a control group. It clearly highlights the potential of fetal ECG as a potential clinical screening tool to aid diagnosis and management of compromised fetuses.

#### FETAL HEART RATES AND VALVE TIMING INTERVALS

Fetal heart rate (FHR) monitoring by Doppler based Cardiotocography (CTG) in the third trimester is a commonly established method to identify fetal compromises (Sandmire and DeMott, 1998). However, sometimes abnormal variability in FHR may not necessarily represent the fetus in distress (Murphy et al., 1991; Vincent et al., 1991). Doppler Ultrasound can provide more information on fetal wellbeing and development beyond the FHR. Using FECG as a reference, automated techniques were proposed in order to identify fetal cardiac valve motion from Doppler signals, providing new measures of mechanical and electromechanical activity of the fetal heart (Marzbanrad et al., 2014, 2016). It provides systolic time intervals (STI) of the fetal cardiac cycle, which have been analyzed by several authors in the past and showed differentiation of fetuses with a variety of perinatal problems (Organ et al., 1980; Koga et al., 2001; Khandoker et al., 2009).

Al-Angari et al. presented a hybrid Empirical Mode Decomposition-Kurtosis method to estimate beat-to-beat fetal heart rate from continuous Doppler signals and then compare with the same from abdominal lead fetal ECG signals.

Jezewski et al. then showed the evidence of equivalence of those two methods (Doppler signals and fetal ECG signals) in terms of recognition of classical FHR patterns such as baseline, accelerations/decelerations, short- and long-term variabilities. These findings might be very useful in clinical settings, as Doppler based FHR monitors are commonly used in routine obstetrics check-ups.

Marzbanrad et al. explained the method on how to estimate fetal cardiac valves' timing intervals from Doppler signals and how those parameters are correlated with fetal gestational development. The proposed method could provide new measures for fetal physiological development.

## ANIMAL MODEL IN PERINATAL MONITORING RESEARCH

Funamoto et al. demonstrated how a fetal mouse model could be utilized to investigate the fetal brain hemorrhage in an ischemia/reperfusion model by using ultrasound B-mode imaging. The use of fetal mice is a novel model for future perinatal research in a variety of fetal stress and compromise.

# CONCLUDING REMARKS

Effective prediction and prevention of fetal stress, compromise or anomalies have now become an emerging research priority. Traditional methods of screening fetal well-being are still popular within the clinical community, however, they demonstrate clear limitations. This research topic briefly presents an overview of original and relevant contributions covering the areas of emerging new signal processing techniques and algorithms enabling early diagnosis of fetal compromises. It is hoped that the proposed technologies and systems could result in improved fetal health management and treatment at the point of need, reduced unnecessary C-section, and the associated economic burden, offering a better quality of early start of life in this world.

# AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

# ACKNOWLEDGMENTS

The guest editors would like to thank all the contributors of the research topic and all reviewers for their thoughtful and valuable comments.

# REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Khandoker, Marzbanrad and Kimura. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Commentary: Computerised interpretation of fetal heart rate during labour (INFANT): a randomised controlled trial

Martin G. Frasch<sup>1</sup> \*, Geraldine B. Boylan2, 3, Hau-tieng Wu<sup>4</sup> and Declan Devane2, 5

<sup>1</sup> Department of Obstetrics and Gynecology, University of Washington, Seattle, WA, United States, <sup>2</sup> Irish Centre for Fetal and Neonatal Translational Research, Cork, Ireland, <sup>3</sup> Department of Paediatrics and Child Health, University College Cork, Cork, Ireland, <sup>4</sup> Department of Mathematics and Department of Statistical Science, Duke University, Durham, NC, United States, <sup>5</sup> School of Nursing and Midwifery, NUI Galway, Galway, Ireland

Keywords: labor, fHRV, EEG, acidemia, sampling rate, FHR, neonate

#### **A commentary on**

Edited by:

Ahsan H. Khandoker, Khalifa University, United Arab Emirates

#### Reviewed by:

Roland H. Hentschel, Universitätsklinikum Freiburg, Germany Stefan Gebhardt, Stellenbosch University, South Africa Faezeh Marzbanrad, Monash University, Australia

> \*Correspondence: Martin G. Frasch mfrasch@uw.edu

#### Specialty section:

This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology

Received: 07 July 2017 Accepted: 06 September 2017 Published: 28 September 2017

#### Citation:

Frasch MG, Boylan GB, Wu H and Devane D (2017) Commentary: Computerised interpretation of fetal heart rate during labour (INFANT): a randomised controlled trial. Front. Physiol. 8:721. doi: 10.3389/fphys.2017.00721

**Computerised interpretation of fetal heart rate during labour (INFANT): a randomised controlled trial**

by The INFANT Collaborative Group (2017). Lancet 389, 1719–1729. doi: 10.1016/S0140-6736(17) 30568-8

The recent study by the INFANT group reported no evidence of benefit on neonatal outcomes associated with the use of decision-support software in conjunction with cardiotocography (CTG) compared with CTG alone (The INFANT Collaborative Group, 2017). Concerns about the study design have been voiced (Keith, 2017).

Earlier studies comparing continuous CTG with intermittent auscultation during labor have shown CTG during labor to be associated with reduced rates of neonatal seizures, largely of unknown long term consequence, but no clear differences in other measures of neonatal mortality and morbidity. However, continuous CTG is associated with increased operative delivery rates (Alfirevic et al., 2017). We now know that supplementing CTG with decision-support software is unlikely to improve outcomes. Yet, given the absence of any alternative to CTG, for high risk women, it is likely that conventional CTG will retain a firm grip in labor wards.

Efforts should now be directed toward interrogating fetal heart rate variability (fHRV) more deeply and/or measuring and evaluating different physiological parameters of intrapartum fetal well-being. For example, the true predictive ability of fetal ECG can only be determined once it is collected at a sampling rate that preserves more of the underlying physiological information. We have demonstrated this approach to assess the severity of hypoxic ischaemic encephalopathy (HIE) in the immediate newborn period (Goulding et al., 2017). Animal and human studies show that the current mode of ECG acquisition is outdated, imprecise, and discards important predictive information (Durosier et al., 2014; Li et al., 2015). To address this challenge, we recently developed and validated an algorithm for low-cost, portable high quality maternal, and fetal ECG monitoring capable of working with one or two maternal abdominal ECG channels to extract the fetal ECG (Li and Wu, 2017; Wu et al., 2017). While we need at least two channels for the algorithm introduced in (Wu et al., 2017), it could be applied to handle the single channel maternal abdominal ECG signal as generalized in (Li and Wu, 2017). An important challenge is refining the algorithm to perform well in the case of twin pregnancies. Also, for the wide acceptance of the technology, it will have to be considerably less expensive than the currently available fetal ECG monitors.

Surprisingly little research has been done on what is perhaps the most obvious and direct source of predictive information about fetal brain health, the Electroencephalogram (EEG), even though it is very sensitive to hypoxia-ischaemia (Murray et al., 2016; Finn et al., 2017). Using animal models, we have developed and validated a new algorithm for fetal EEG as predictor of academia (Frasch et al., 2015) and are beginning a clinical study using a fetal EEG prototype device (NCT03013569). Our data shows that pathognomonic EEG changes emerge at a pH ≤ 7.20, which the INFANT study reported in at least 12% of all babies.

We expect that direct acquisition of fetal EEG during labor and a more precise acquisition of fetal ECG (and fHRV) offers potential in preventing acidosis and brain injury due to intrapartum hypoxia-ischaemia. Fetal infection is an important contributor to perinatal brain injury and fHRV has been shown to reflect fetal inflammatory response systematically and in an organ-specific manner (Durosier et al., 2015; Frasch et al., 2016; Liu et al., 2016). Prospective studies in large anteand intrapartum cohorts of pregnant women and fetuses are now needed to validate the utility of fetal ECG and EEG

#### REFERENCES


monitors. Such studies will not only yield biomarkers to predict complications during delivery due to underlying infection, hypoxia or acidemia, but also pinpoint the risks for abnormal postnatal developmental trajectories. Knowing on which babies to focus will provide a foundation upon which to build the therapeutic strategies to correct early deviations from healthy developmental trajectories.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

MF acknowledges support by CIHR and FRQS. HW is partially supported by Sloan Research Fellow FR-2015-65363.


**Conflict of Interest Statement:** MF is an inventor of related patent application entitled "EEG Monitor of Fetal Health" including U.S. Patent 9,215,999. MF and HW filed a patent for the aECG method referred to in this commentary.

The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Frasch, Boylan, Wu and Devane. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Fetal Cardiac Doppler Signal Processing Techniques: Challenges and Future Research Directions

*Saeed Abdulrahman Alnuaimi1 \*, Shihab Jimaa1 and Ahsan H. Khandoker <sup>2</sup>*

*1Department of Electrical and Computer Engineering, Khalifa University, Abu Dhabi, United Arab Emirates, 2Department of Biomedical Engineering, Khalifa University, Abu Dhabi, United Arab Emirates*

#### *Edited by:*

*Joseph L. Greenstein, Johns Hopkins University, United States*

#### *Reviewed by:*

*Jason H. Yang, Massachusetts Institute of Technology, United States Clarissa Lim Velayo, University of the Philippines Manila, Philippines Radek Martinek, Technical University of Ostrava, Czechia*

#### *\*Correspondence:*

*Saeed Abdulrahman Alnuaimi saeed.alnuaimi@kustar.ac.ae*

#### *Specialty section:*

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Bioengineering and Biotechnology*

*Received: 21 May 2017 Accepted: 11 December 2017 Published: 22 December 2017*

#### *Citation:*

*Alnuaimi SA, Jimaa S and Khandoker AH (2017) Fetal Cardiac Doppler Signal Processing Techniques: Challenges and Future Research Directions. Front. Bioeng. Biotechnol. 5:82. doi: 10.3389/fbioe.2017.00082*

The fetal Doppler Ultrasound (DUS) is commonly used for monitoring fetal heart rate and can also be used for identifying the event timings of fetal cardiac valve motions. In early-stage fetuses, the detected Doppler signal suffers from noise and signal loss due to the fetal movements and changing fetal location during the measurement procedure. The fetal cardiac intervals, which can be estimated by measuring the fetal cardiac event timings, are the most important markers of fetal development and well-being. To advance DUS-based fetal monitoring methods, several powerful and well-advanced signal processing and machine learning methods have recently been developed. This review provides an overview of the existing techniques used in fetal cardiac activity monitoring and a comprehensive survey on fetal cardiac Doppler signal processing frameworks. The review is structured with a focus on their shortcomings and advantages, which helps in understanding fetal Doppler cardiogram signal processing methods and the related Doppler signal analysis procedures by providing valuable clinical information. Finally, a set of recommendations are suggested for future research directions and the use of fetal cardiac Doppler signal analysis, processing, and modeling to address the underlying challenges.

Keywords: fetal Doppler, signal processing, fetal cardiac intervals, fetal monitoring, fetal heart rate, fetal cardiology

# INTRODUCTION

Fetal heart rate (FHR) monitoring has been extensively used to assess fetal well-being. The process of FHR monitoring is commonly used during prenatal screening to detect possible fetal health problems that may result in neurological damage or in some cases fetal death during labor. Statistics have shown that 1 out of every 125 babies is born with some kind of congenital cardiac defect (Anisha et al., 2014). When certain pregnancy risk factors have been identified, the FHR must be monitored during labor as a routine physiological measurement (Elmansouri et al., 2014). Cardiotocography (CTG) is the standard methodology in hospitals to monitor fetal well-being and is based on the recording of FHR using a special device. This methodology is popular and commonly used because it is easy to use, is a non-invasive technique, has no contraindications, and can be used frequently. The main disadvantage of CTG is its high sensitivity to fetal movement, as the detection of FHR mostly relies on the correct positioning of the ultrasound probe. This probe therefore needs to be adjusted in the case of fetal movement to afford accurate measurements. Furthermore, the Doppler ultrasound (DUS) transducer is uncomfortable, and the FHR monitoring procedure involves sending a 2 MHz signal toward the fetus. Consequently, it is considered an invasive method and it is not recommended, especially for recordings over long periods under severe conditions (Maeda, 1990; Elmansouri et al., 2014).

In this review paper, the fetal electrocardiogram (FECG) is mentioned many times as a reference to compare and verify the results of the latest fetal cardiac Doppler signal processing techniques. The FECG carries essential information about the fetal heart function. The characteristics or features of the FECG are vital for revealing the fetal development, as well as the existence of fetal distress or congenital heart defects. The FECG is considered to be an effective tool for diagnosing specific structural defects. Currently, there are two methods for recording the FECG: indirect and direct measurements. ST analysis is considered to be a direct method and is performed by directly attaching an electrode to the scalp of the fetus to provide a clean ECG signal. However, this method is not used extensively because of its inherent danger to both the mother and fetus and it can be used in clinical practice only after 36 weeks of gestation. Indirect measurement involves non-invasive FECG recordings that are obtained by placing a skin electrode on the mother's abdomen. Fetal monitoring is based entirely on the FHR because there is no available technology to reliably measure FECG (Chandraharan, 2010; Anisha et al., 2014; Bsoul, 2015).

Most fetal cardiac defects have some variation in their morphology, which reflects the health status of the fetus, although morphological analysis is lacking in other conventional methods. Most of the clinically essential data in the FECG signal are embedded in the amplitude and duration of its waveforms. Fetal cardiac waveforms help doctors to diagnose fetal arrhythmias such as bradycardia, tachycardia, asphyxia, and congenital heart disease (Voicu et al., 2010; Anisha et al., 2014).

This paper is organized as follows: Section "Introduction" discusses the principles of Doppler-based fetal monitoring, the challenges of fetal cardiac Doppler monitoring, and its biomedical applications. Section "Fetal Heart Rate Monitoring" provides an overview about FHR monitoring and a comprehensive survey of the fetal cardiac Doppler methods based on signal processing techniques. Section "Future Work" provides a summary of the existing challenges and suggests potential directions for future research.

#### Fetal Cardiac Doppler Ultrasound

Although fetal DUS was discovered many years ago, research interest in this field has only arisen over the last few years. **Figure 1** shows the number of articles published in the Institute of Electrical and Electronics Engineers (IEEE), ScienceDirect, and the National Institutes of Health (PubMed) databases. The search was focused on both fetal cardiac Doppler and CTG specifically. **Figure 2** shows a simultaneously captured fetal cardiac DUS and FECG. The high-frequency component of the DUS signal allows the manual identification of the fetal cardiac valve movements. These sensitive markers can be used for assessing the fetal cardiac performance (Merz, 2004; Marzbanrad et al., 2014).

The form of the fetal cardiac Doppler signal changes over time, as a result of the changes in the transducer location in relation to the moving signal source. This non-stationary nature results in changes in the fetal cardiac Doppler signal on a beatto-beat basis, in addition to alteration of the spectral characteristics and signal content over time. The resulting complexity makes precise measurement of time intervals in the cardiac cycle quite challenging. The accuracy of calculating the FHR using DUS is very low, and the manual identification of beat-tobeat opening and closing of valves is time consuming, requires special expertise, and is subject to inter- and intra-observer and visual errors. This challenge can be solved by the automated identification of valve movement. Furthermore, DUS signal quality is crucial for the reliable estimation of the fetal cardiac valve timings and also to provide real-time feedback to the operator during data collection; this issue can be solved using automated DUS quality assessment (Wrobel, 2004; Marzbanrad et al., 2015a, 2016).

# Biomedical Applications of Doppler Ultrasound

#### Echocardiography

An echocardiogram, often known as a cardiac echo, is a sonogram of the heart. It allows the visualization of the internal structure of the fetal heart by creating images using the standard 2D, 3D, and Doppler ultrasound. This technique is used routinely during pregnancy in the diagnosis, management, and follow-up of pregnant women carrying a fetus with any suspected or known heart abnormality. This advanced method can provide a significant amount of helpful information, such as pumping capacity, the shape of the heart, and the presence and location of any cardiac tissue damage. Furthermore, it can also allow physicians and researchers to perform a wealth of helpful calculations to estimate the heart function, including the calculation of the cardiac output and heart diastolic function, which indicates how well the heart relaxes.

The most important advantage of echocardiography is that it is a non-invasive technique without any major side effects, except being exposed many times to its high-frequency sound waves which might harm the patient. Pulsed-wave Doppler ultrasound can be used to generate an accurate assessment of the blood flowing through the heart. To visualize any abnormalities, color Doppler and spectral Doppler can be used. Moreover, these methods can be used for discovering any valve-related issues, such as blood leaking through the valves and estimating how well the valves open (Kwon and Park, 2016; Michelfelder et al., 2016).

#### Doppler Fetal Monitor

In 1964, Callagan developed the Doppler ultrasound monitor. A Doppler FHR monitor is a handheld ultrasound transducer used to detect the fetal heart beat during prenatal care. This device uses the Doppler effect to provide an audible simulation of the heartbeat. The use of this monitor is sometimes known as Doppler auscultation. To enhance the sonography, Doppler measurements can be made by employing the Doppler effect to assess the movement of the cardiac structures toward or away from the probe and the corresponding velocity of this movement. For example, the speed and direction of the blood flow in a fetal heart valve can be determined and visualized by calculating the frequency shift of a particular sample volume. The current advanced ultrasound transducers use pulsed Doppler ultrasound for the velocity measurements. The pulsed-wave scanners operate by transmitting and receiving a series of pulses, and the frequency shift can be obtained using the relative phase changes of the pulses. The most important advantage of pulsed Doppler over continuous-wave Doppler is that it permits distance information to be obtained and gain correction to be applied. However, the disadvantage is that the Doppler measurements can be affected by aliasing (Feinstein et al., 1993; Jezewski et al., 2008).

# FHR MONITORING

Fetal heart rate monitoring provides important information regarding fetal cardiac conditions and is a standard clinical procedure that is widely used during pregnancy and labor to assess the well-being of the fetus. FHR monitoring has become a common clinical practice since its introduction at Yale University in 1958. FHR tracing can be performed automatically using electrocardiograms and CTG, whereas FHR monitoring and its interpretation are conducted manually by obstetricians, which can lead to substantial inter-observer and intra-observer differences. In order to allow accurate and consistent clinical practice, automated FHR monitoring and interpretation is recommended (Hon, 1958).

A number of methods have been developed to address various issues in computerized automatic FHR monitoring and interpretation, including FHR signal modeling and representation, feature extraction, and pattern classification. In FHR monitoring devices, the Doppler effect is the preferred technique. Furthermore, a reliable fetal heart simulator becomes essential for testing Doppler FHR monitoring devices (Liu et al., 2008).

## Developing a Doppler FHR Testing Device and Generating the Doppler Frequency Shift

A number of authors have designed a device that can simulate the fetal heart valves and cardiac wall motion in air. This device can be used to test Doppler FHR monitoring in a clinical environment by using a modified electrical relay in air and generating a similar Doppler frequency shift to the fetal heart activity. Accordingly, by modifying the opening and closing velocities of the relay, similar Doppler frequency shifts to the fetal heart activity in soft tissue detected by the ultrasound probe can be produced. The authors tested the functionality and accuracy of the new device with current commercially available Doppler FHR monitors by focusing on the effect of the relay on the frequency range of the Doppler shift. A band-pass filter was used to eliminate the effects of fetal breathing movement, hiccup movement, global movement, and maternal activity, to allow calculation of the periodicity of the frequency reflected by fetal cardiac activity. The authors provided a comparison between the Doppler frequency shifts of this device and the measurements from previous studies and its applied filters to describe the Doppler frequency shift by fetal cardiac activity given in Table S3 in Supplementary Material (Mert et al., 2015).

After comparing the estimated Doppler frequency shifts of the new device with the previously proposed techniques, both the Doppler shift of the fetal cardiac valve and the wall in the same period could be simulated by determining the difference between the opening and closing velocities of the relay armature. The proposed method was tested with commercially available Doppler FHR monitors and the authors claimed that only 4 tests out of 10 displayed a trivial error. Therefore, the proposed device functionality and precision have been proved, demonstrating the potential accuracy of Doppler FHR monitoring devices and the possibility of avoiding any errors in the clinical diagnosis procedures. The authors also claimed that the easy and low-cost implementation could make this device a good candidate for further applications in the future (Mert et al., 2015).

# Real-Time Signal Processing for FHR Monitoring

A new algorithm has been developed based on adaptive thresholding, differencing of local maxima and minima, digital filtering, and statistical properties in the time domain. This algorithm can be used to simultaneously measure the FHR and maternal heart rate from the maternal abdominal electrocardiogram during both pregnancy and labor for ambulatory monitoring. The researchers used a microcontroller-based system to implement the proposed algorithm in the real-time domain. For statistical comparison, a Doppler ultrasound fetal monitor was used on five volunteers with low-risk pregnancies at a gestational age of 35–40 weeks. The proposed algorithm was reported to deliver an average percent root-mean-square difference of 5.32% and a linear correlation coefficient of 0.84–0.93 (Ibrahimy et al., 2003). Furthermore, the FHR curves remained inside a ±5 beats per minute limit relative to the reference ultrasound method for 84.1% of the time (Ibrahimy et al., 2003).

### A Novel Technique for FHR Estimation from the DUS Signal

A new method was recently proposed based on providing accurate beat-to-beat cardiac cycle values of the FHR through multiple measurements of a given fetal cardiac cycle in the DUS signal. The proposed algorithm involves three stages: (i) the dynamic adjustment of the autocorrelation window, (ii) the adaptive autocorrelation Doppler signal peak detection, and (iii) the determination of beat-to-beat time durations. The researchers compared the estimated FHR values and the calculated indices describing the variability of FHR to the reference data obtained from the direct fetal ECG, as well as to another method for FHR estimation. The authors claimed that the results showed that the proposed method increases the estimation precision in comparison to the conventionally used FHR monitoring devices, and this enabled them to calculate reliable parameters describing the FHR variability. Furthermore, according to a comparison of the results of the proposed method to those of the other FHR estimation methods applied, the former technique rejected a much lower number of measured cardiac cycles as being unacceptable (Jezewski et al., 2011).

The proposed technique for the FHR estimation on a beat-tobeat basis provided a high accuracy for the fetal heart interval measurement and enabled a reliable, accurate, and quantitative assessment of the FHR variability, while reducing the number of invalid cardiac cycle measurements. Comparing the proposed method with Peters' method (CHL et al., 2004), the authors noted a very similar accuracy for the measured intervals and slightly superior findings in terms of the variability indices of the single-measurement method (mean relative error: −5.1%). The authors also mentioned the requirement for assessing the short-term variability for reliable evaluation of FHR signals because it does not depend directly on the accuracy of the fetal cardiac interval measurement. However, in the case of both the proposed method and Peters' method, the mean relative errors of the short-term variability did not exceed −7%, which might be a satisfactory result (Jezewski et al., 2011).

## Fetal Cardiac Systolic Time Intervals (STIs)

Several antenatal fetal cardiac assessment techniques have been developed to evaluate antepartum fetal cardiac risks. These represent sensitive markers for assessing fetal cardiac performance and permit the evaluation of the fetal cardiac electromechanical coupling. This evaluation process is a fundamental and clinically significant aspect of determining the heart physiology (Weissler et al., 1968; Lewis et al., 1977; Marzbanrad et al., 2014, 2015a).

The main basis for estimating these electromechanical markers are the opening and closing timings of the fetal cardiac valves. **Figure 3** shows the STIs, which are one of the markers that have received considerable attention as an indicator of myocardial function. From a clinical perspective, the most convenient of the STIs are the pre-ejection period (PEP), the isovolumetric contraction time (ICT), and the left ventricular ejection time. The PEP is a sensitive sign of the functional state of the fetal myocardium and becomes prolonged early in the development

of hypoxemia and acidosis. The ICT has also been suggested as a reliable index to represent fetal cardiac contractility. To obtain and analyze the STIs, several non-invasive methods have been proposed such as DUS and abdominal ECG. The cardiac timings were identified manually after filtering the DUS signal using a band-pass filter. The poor quality of the abdominal ECG and the high variability of the fetal cardiac Doppler signal over time are the main limitations of this technique (Murata and Martin, 1974; Murata et al., 1978; Koga et al., 2001; Shakespeare et al., 2001; Yumoto et al., 2005; Marzbanrad et al., 2014).

The DUS signal was divided into different frequency shift ranges by using a digital narrow band-pass filter. From the peaks in one of the filtered signals, the mitral and aortic valve motions were identified. In contrast, other researchers used the short-time Fourier transform method (STFT) to analyze the DUS signal. Correspondingly, the high-frequency component of the DUS signal is related to the valve movements, whereas the low-frequency component is linked with the cardiac wall motion (Koga et al., 2001; Shakespeare et al., 2001).

## The DUS Signal Multiresolution Wavelet Analysis

In another study (Khandoker et al., 2009a), the multiresolution wavelet analysis technique was applied to the fetal cardiac Doppler signal (Marzbanrad et al., 2016). Wavelet analysis is a powerful technique to decompose non-stationary signals with variable spectral characteristics over time. Using this technique, the fetal cardiac DUS signal is first decomposed into different scales with corresponding resolution levels. From the results, it is clear that the multiresolution wavelet analysis enabled the frequency contents of the fetal cardiac Doppler signals to be linked to the opening (o) and closing (c) of the fetal heart valves [aortic (A) and mitral (M)]. This technique was tested on DUS signal samples at a gestational age of 28–36 weeks. The valve movements were visualized as peaks in the detailed DUS signal at level 2 wavelet decomposition. The next stage is to assign each peak manually to the opening and closing of the cardiac valves.

#### Blind Source Separation (BSS) Method

The BSS method was used to extract the fetal ECG from the abdominal ECG mixture, because the abdominal ECG is noisy and observing the fetal cardiac R wave is very difficult. Moreover, the correlation of the fetal cardiac cycle R–R interval with the interval of the R wave to each valve motion was examined, which has potential clinical applications. This correlation was introduced as a criterion for diagnosing fetal cardiac abnormalities. Previous researchers have investigated the automatic identification of these abnormalities as stated above (Sato et al., 2007; Khandoker et al., 2009b, 2010; Marzbanrad et al., 2014).

A non-invasive and automated method using an integrated fetal transabdominal ECG system and Doppler cardiogram was developed for the identification of fetal cardiac abnormalities (Khandoker et al., 2010). The authors used both the multiresolution wavelet analysis method and the Jensen–Shannon divergence method for the identification of the frequency contents of the Doppler signals for subsequent linking to the opening and closing of the fetal heart valves (aortic and mitral). Until recently, these cardiac time intervals have mostly been measured by ultrasound with manual identification of the fetal cardiac valve motion events. The potential for automated non-invasive assessment without obstetric ultrasonography expertise allows the development of integrated fetal ECG and cardiac Doppler signals from heart valve and wall motion events. Various applications of this technology may be feasible, enabling assessment of its value in antenatal monitoring of the fetal heart well-being (Khandoker et al., 2010; Marzbanrad et al., 2015a).

## Automated Estimation of Fetal Cardiac Timing Events from the DUS Signal Using the Hybrid Support Vector Machine (SVM)/ Hidden Markov Model (HMM) Model

An automated methodology has been proposed to identify the occurrence of the cardiac events based on the sequence of the movements, timings, and patterns of both the valve and the wall in the fetal DUS signal components. These researchers proposed using the empirical mode decomposition (EMD) instead of STFT or wavelet analysis. This technique has been widely used for many applications, such as image processing, speech processing, and biomedical signal processing. The authors introduced three approaches to be combined with the EMD method for automated identification: the HMM, SVM, and hybrid SVM/ HMM. In this case, the changes of the cardiac intervals were evaluated from the 16th to the 41st week of gestation (Echeverria et al., 2010; Mijovic et al., 2010; Marzbanrad et al., 2014, 2015a; Springer et al., 2016).

By comparing the results by Marzbanrad et al. (2014) with the pulsed Doppler image by Khandoker et al. (2009a), it can be concluded that the proposed hybrid algorithm afforded better results in terms of identification of the cardiac events. Moreover, compared to the conventional manual identification method by Khandoker et al. (2009a), a higher percentage of the fetal cardiac valve movement events was identified. In the study by Marzbanrad et al. (2014), the FECG was used as a reference for segmentation to facilitate the estimation of the timing of cardiac events. The proposed algorithm results provide a continuous and beat-to-beat identification of the fetal cardiac intervals, which can be used later for clinical purposes. Moreover, the authors noted a limitation in the proposed algorithm, in that a quantitative comparison with the pulsed-wave Doppler imagebased valve motion timings was not provided. The more accurate method of trans-vaginal pulsed Doppler imaging can be used for fetuses in the first trimester. However, the authors claimed that the proposed method is compatible with wide continuous FHR monitoring during the second to third trimesters. The authors also stated that the need for more precise quantitative comparison of the proposed method results with pulsed Doppler images will require advanced research, such as image processing and recognition processes, which are beyond the scope of this study. Moreover, the authors suggested using the quantitative comparison in future studies (Kikallio et al., 2005; Marzbanrad et al., 2014).

# Identifying the Cardiac Events Automatically Using Machine Learning Techniques

In the conventional technique for identifying fetal cardiac events, the fetal cardiac timing events were manually assigned to the signal peaks and the time intervals were calculated. According to a previous study (Marzbanrad et al., 2014), the aim is to identify the cardiac events automatically using machine learning techniques. In this regard, each peak in the DUS signal should be classified as an indicator of one of the fetal cardiac valve timing events or none of them. The aforementioned researchers (Marzbanrad et al., 2014) used the HMM; this statistical model can allow the fetal cardiac events to be determined based on the probabilistic model of their occurrence sequence and timings in the DUS signal. However, the authors claimed that the amplitude as well as the timing of the peaks can also be used for the event classification.

Moreover, the authors used a powerful classifier, namely, SVM, for the classification of the fetal cardiac events. The authors noted certain limitations of SVM, including not considering the temporal dependence of the occurrence of events, faults in peak classification in some cardiac cycles, or an incorrect order of cardiac events. Furthermore, they proposed the hybrid HMM-SVM method as a solution to overcome the limitations of SVM and HMM. Each fetal cardiac cycle was segmented with reference to the FECG. The results of the proposed method showed the following achievements: 94.5% of mitral opening, 91.1% of mitral closing, 95.3% of aortic valve opening, and 98.8% of aortic valve closing. The authors analyzed the changes of the fetal cardiac intervals for three fetal age groups: 16–29, 30–35, and 36–41 weeks. These timings were identified using the proposed hybrid HMM-SVM technique, which provided more accurate results than the conventional manual approaches. The pulsed Doppler images were then used to verify the identified timings.

The K-means clustering method was also used to find the fetal DUS cardiac component patterns and match the DUS components of each cardiac cycle beat-to-beat to one of the models. To decompose the non-stationary signals with variable spectral characteristics over time, the authors proposed a multiresolution wavelet analysis method to apply to the DUS signal. The fetal cardiac valve movements were visualized as peaks in the detailed signal by using the wavelet analysis at level 2 wavelet decomposition. Subsequently, the hybrid SVM-HMM was used to identify the fetal cardiac valve motion events from the peaks of the DUS component and the model was trained specifically for its corresponding cluster (Marzbanrad et al., 2016).

# Hybrid EMD-Kurtosis Method

A new method has been proposed to estimate FHR and its variability from fetal DUS based on EMD and kurtosis statistics (Al-Angari et al., 2017). This method relies on computing the kurtosis function on the intrinsic mode functions extracted from the DUS signal to estimate cardiac beat-to-beat intervals. The authors also provided a comparison between the estimated beat-to-beat intervals using the proposed method and the autocorrelation function (AF) with respect to the R–R intervals computed from FECG. This method was tested on DUS signals from 44 pregnant mothers in the early (20 cases) and late (24 cases) gestational weeks. The authors reported that the EMD-kurtosis method showed superior performance for estimating the mean beat-tobeat intervals, with an average difference of 1.6 ms from the true mean R–R intervals, compared with the value of 19.3 ms obtained using the AF method. The EMD-kurtosis method is more robust than AF for low SNR cases and can be used in a hybrid system to estimate the beat-to-beat intervals from the fetal DUS signal. A limitation of this method is that the collected data were 1 min in length, which is insufficient to control the fetal states. This method can be easily implemented in a microprocessor in the same DUS machine or in a separate device for practical clinical use (The Society for Research in Child Development, 2015; Apostolidis and Hadjileontiadis, 2016).

In this study, the authors detected six different patterns for the DUS component. Two major findings were obtained. First, the generated patterns 1 and 6 occurred with significantly higher rates for the age group after 36 weeks than for the age group before 32 weeks. Second, the remaining patterns 3, 4, and 5 were observed with significantly higher rates for the early gestation group, whereas the percentage difference between the two age groups was not significant for pattern 2. The authors claimed that by comparing the proposed method with the technique in which clustering was not performed, the former allowed better precision and recall to be achieved. It was reported that the occurrence rates of five of the cardiac patterns differed between the fetuses older than 36 weeks and those under 32 weeks. Moreover, each cardiac pattern exhibited its own characteristics, such as amplitude range and timing of the peaks linked to the aortic and mitral valve motion. The authors compared the clustering method results to improve the identification of opening and closing of the mitral valve by the SVM-HMM method to the method without clustering, as verified using pulsed-wave Doppler images. The average precision and recall of the proposed technique with clustering were higher than those of the method without clustering by 83.4 and 84.2%, respectively (Marzbanrad et al., 2016).

In this review paper, several techniques have been discussed to address the challenges in fetal cardiac Doppler signal processing. These techniques facilitate the DUS signal analysis and classification, linear decomposition, and automated identification of the fetal cardiac events using machine learning techniques. Multiresolution wavelet analysis and BSS are techniques used to decompose non-stationary signals with variable spectral characteristics over time. This analysis has enabled the frequency contents of the fetal cardiac Doppler signals to be linked to the opening (o) and closing (c) of the fetal heart valves (aortic and mitral). These techniques have not yet been automated in terms of assigning each peak to the opening and closing of the cardiac valves.

Conventionally, fetal cardiac timing events are manually assigned to the signal peaks and the time intervals are calculated. This process requires experts and consumes time. The latest approach in fetal Doppler signal processing is the automated identification of the occurrence of the cardiac events based on the sequence of the movements, timings, and patterns of both the valve and the wall in the fetal DUS signal components. This technique consists of three approaches to be combined with the EMD method for automated identification: the HMM, SVM, and the hybrid SVM/HMM method. The technique was used to analyze the changes in the fetal cardiac intervals for three fetal age groups: 16–29, 30–35, and 36–41 weeks. Although this technique led to an improvement in the estimation of fetal cardiac timing events, more precise quantitative comparison of the results is required. In addition, the K-means clustering method has been used to ascertain the fetal DUS cardiac component patterns and match the DUS components of each cardiac cycle "beat-to-beat" to one of the models. Six different patterns for the DUS component were found. The use of the clustering method improved the identification of the opening and closing of the mitral valve by the SVM-HMM method compared with the method without clustering. The key limitation of the current automated estimation of fetal cardiac intervals from DUS signals is the complexity of these systems.

# FUTURE WORK

Several issues concerning fetal cardiac Doppler signal processing still need to be resolved. First, one limitation of the latest new technique for computerized estimation of fetal cardiac time intervals from DUS signals is that the quantitative comparison with the pulsed-wave Doppler image-based valve motion timings has not yet been achieved (Marzbanrad et al., 2015a). Trans-vaginal pulsed Doppler imaging is considered one of the most accurate techniques, yet it can only be performed in the first trimester. The latest computerized estimation of fetal cardiac time intervals from the fetal cardiac DUS signal will help in continuing with the same accuracy during the remainder of the pregnancy. As for future studies, previous researchers suggested the application of image recognition processes on the pulsed Doppler images to enhance the quantitative comparison of the result accuracy (Marzbanrad et al., 2014).

In order to measure the correctness of the automated identification of the fetal cardiac valve motion by hybrid SVM-HMM methods, M-mode and pulsed-wave Doppler images were used (Marzbanrad et al., 2014, 2016). The findings suggested that more quantitative comparisons should be performed in future investigations. As the fetus grows in the late gestation period, its movement increases, which results in the non-stationary nature of the fetal cardiac Doppler signal. Longer fetal cardiac Doppler recordings are required to evaluate the non-stationary nature and improve the FHR accuracy. The advanced analysis of the fetal cardiac non-stationary Doppler signal and its influence on the FHR monitoring remains a topic for future studies (Marzbanrad et al., 2015b). The variation of cardiac features is a limitation for the predefined range of the valve motions in the

#### REFERENCES


proposed method (Marzbanrad et al., 2015a). The recommendation of these researchers for future studies is to assess the validity of the measurements on the abnormal cases. Finally, reducing the complexity of the current automated estimation of fetal cardiac intervals from DUS signals should be subject to future investigation.

# CONCLUSION

The current survey has discussed the importance of clinical FHR monitoring and provided a critical overview of the existing techniques. Several solutions have been proposed to overcome the effect of the non-stationary nature of the fetal cardiac Doppler signal on a beat-to-beat basis, in addition to the variation of the spectral characteristics and signal content over time. These solutions include band-pass filters, STFT, and wavelet analysis. With respect to DUS signal analysis and classification, linear decomposition techniques such as BSS for source separation have shown promising results in the isolation of fetal cardiac Doppler signal components. Furthermore, several automated techniques for the identification of valve movements using fetal cardiac Doppler signals have been discussed to overcome the shortcomings of manual techniques, including their time-consuming nature, and to increase the FHR measurement accuracy. It is important to mention that using machine learning techniques within the context of FHR monitoring can potentially open up new research areas that have so far gone unexplored. Furthermore, potential directions for future research have been suggested to improve current fetal cardiac Doppler signal processing and analysis techniques. Finally, the growing interest in fetal cardiac DUS is set to provide new opportunities for reliable and accurate FHR monitoring.

# AUTHOR CONTRIBUTIONS

SA is the first author who is a PhD student, AK is the main supervisor, and SJ is the secondary supervisor.

### FUNDING

This study was supported by research incentive funds from Khalifa University Internal Research Fund Level 2.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://www.frontiersin.org/articles/10.3389/fbioe.2017.00082/ full#supplementary-material.


closing timings in developing human fetuses. *IEEE J. Biomed. Health Inform.* 20, 240–248. doi:10.1109/JBHI.2014.2363452


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2017 Alnuaimi, Jimaa and Khandoker. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Template-based Quality Assessment of the Doppler Ultrasound Signal for Fetal Monitoring

Camilo E. Valderrama<sup>1</sup> \*, Faezeh Marzbanrad<sup>2</sup> \*, Lisa Stroux <sup>3</sup> and Gari D. Clifford4, 5

<sup>1</sup> Department of Mathematics and Computer Science, Emory University, Atlanta, GA, United States, <sup>2</sup> Department of Electrical and Computer Systems Engineering, Monash University, Melbourne, VIC, Australia, <sup>3</sup> Department of Engineering Science, Institute of Biomedical Engineering, University of Oxford, Oxford, United Kingdom, <sup>4</sup> Department of Biomedical Informatics, Emory University, Atlanta, GA, United States, <sup>5</sup> Department of Biomedical Engineering, Georgia Institute of Technology, Atlanta, GA, United States

#### Edited by:

Joseph L. Greenstein, Johns Hopkins University, United States

#### Reviewed by:

Radek Martinek, Technical University of Ostrava, Czechia Paolo Melillo, Università degli Studi della Campania "Luigi Vanvitelli" Caserta, Italy

#### \*Correspondence:

Camilo E. Valderrama cvalder@emory.edu Faezeh Marzbanrad faezeh.marzbanrad@monash.edu

#### Specialty section:

This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology

> Received: 20 April 2017 Accepted: 04 July 2017 Published: 18 July 2017

#### Citation:

Valderrama CE, Marzbanrad F, Stroux L and Clifford GD (2017) Template-based Quality Assessment of the Doppler Ultrasound Signal for Fetal Monitoring. Front. Physiol. 8:511. doi: 10.3389/fphys.2017.00511 One dimensional Doppler Ultrasound (DUS) is a low cost method for fetal auscultation. However, accuracy of any metrics derived from the DUS signals depends on their quality, which relies heavily on operator skills. In low resource settings, where skill levels are sparse, it is important for the device to provide real time signal quality feedback to allow the re-recording of data. Retrospectively, signal quality assessment can help remove low quality recordings when processing large amounts of data. To this end, we proposed a novel template-based method, to assess DUS signal quality. Data used in this study were collected from 17 pregnant women using a low-cost transducer connected to a smart phone. Recordings were split into 1990 segments of 3.75 s duration, and hand labeled for quality by three independent annotators. The proposed template-based method uses Empirical Mode Decomposition (EMD) to allow detection of the fetal heart beats and segmentation into short, time-aligned temporal windows. Templates were derived for each 15 s window of the recordings. The DUS signal quality index (SQI) was calculated by correlating the segments in each window with the corresponding running template using four different pre-processing steps: (i) no additional preprocessing, (ii) linear resampling of each beat, (iii) dynamic time warping (DTW) of each beat and (iv) weighted DTW of each beat. The template-based SQIs were combined with additional features based on sample entropy and power spectral density. To assess the performance of the method, the dataset was split into training and test subsets. The training set was used to obtain the best combination of features for predicting the DUS quality using cross validation, and the test set was used to estimate the classification accuracy using bootstrap resampling. A median out of sample classification accuracy on the test set of 85.8% was found using three features; template-based SQI, sample entropy and the relative power in the 160 to 660 Hz range. The results suggest that the new automated method can reliably assess the DUS quality, thereby helping users to consistently record DUS signals with acceptable quality for fetal monitoring.

Keywords: doppler ultrasound, empirical mode decomposition, fetal monitoring, dynamic time warping, sample entropy, signal quality

# 1. INTRODUCTION

Although medical care has reduced mortality rates across the globe, birth has still remained an event of extreme risk. Approximately 2.6 million stillbirths and 2.8 million early neonatal deaths occur each year (World Health Organization, 2016a,b). Different factors contribute to this high burden, such as the lack of specialized medical professionals and the high cost of the medical devices, mainly affecting low and middle income countries (LMICs). Leading causes for perinatal mortality include Intrauterine Growth Restriction (IUGR) and congenital abnormalities of which, Congenital Heart Disease (CHD) is the most common (van der Linde et al., 2011; Gardosi et al., 2013; Lawn et al., 2016). These abnormalities are currently being detected using ultrasound imaging and more specifically, fetal echocardiography is performed for CHD diagnosis. However, these techniques are expensive and can only be performed by trained sonographers or physicians; hence, their use is limited in LMICs (McClure et al., 2014).

Due to the high incidence and fatal consequences of these abnormalities in low-resource settings, affordable perinatal monitoring solutions are required. One of the most widely used, yet affordable methods for perinatal screening is fetal heart rate (FHR) monitoring. This technique has contributed to reduce perinatal and maternal risks through identification of nonreassuring fetal status (Ayres-de Campos and Bernardes, 2010). Moreover, FHR has the potential for detecting IUGR (Nijhuis et al., 2000; Ferrario et al., 2009), as well as CHD complications (Cullen et al., 1992; Berghella et al., 2001).

FHR monitoring is commonly performed through Cardiotocography (CTG) based on one dimensional Doppler Ultrasound (DUS), that is also used in low cost (under \$17) hand-held devices which can be operated by non-experts. This DUS-based low-cost device has been used to develop affordable perinatal monitoring systems, thus facilitating screening in LMICs. Stroux et al. introduced a mobile-health monitoring system, based on a low-cost transducer and operated by illiterate birth attendants, to detect fetal compromise, such as IUGR, in rural Guatemala (Stroux, 2016; Stroux et al., 2016). DUS can also provide more information beyond the FHR, such as the cardiac valve function, which can further facilitate detection of CHDs and assessment of the fetal development (Marzbanrad, 2015; Marzbanrad et al., 2016b).

Despite the benefits of 1D-DUS, it is susceptible to noise affecting its quality, and it is non-stationary due to the fetal movements, which can complicate the FHR monitoring. Since the quality of the recorded signals is critical to properly detect FHR abnormalities, the assessment of the signal quality is an essential part of the recording process. Stroux and Clifford reported that the accuracy of FHR analysis depends on the signal quality, hence the quality should be ensured during the data collection (Stroux and Clifford, 2013). Magenes et al. also showed the necessity of removing CTG signals with low quality before applying methods for detecting fetal anomalies (Magenes et al., 2001). The quality assessment process, enables providing feedback to the operator during data collection, allowing them to retake or exclude the low-quality signals.

To date, little work has been published concerning the quality assessment of the DUS signals. Stroux and Clifford proposed a method to validate the quality of DUS signals recorded using a hand-held device connected to a smart phone (Stroux and Clifford, 2014; Stroux, 2016). For this purpose, they extracted features based on sample entropy, wavelet decomposition coefficients, and the phone's triaxial accelerometer output. To assess the quality, a logistic regression and a support vector machine (SVM) were trained to classify the recordings into noisy and clean categories. The logistic regression model was able to classify the signal quality with an accuracy of 95.14% on test data, while the SVM achieved an accuracy of 94.44%. Marzbanrad et al. proposed an automated method to assess the DUS signal quality for the application of fetal valve motion identification (Marzbanrad, 2015; Marzbanrad et al., 2015). In their method, DUS signals were segmented into cardiac cycles using noninvasive fetal electrocardiogram (fECG) as reference. Then, 12 features including power, statistical and entropy-based measures, were extracted from a frequency range associated with the fetal cardiac valve motion. Using these features, the signals collected from 57 fetuses were classified as good and poor quality, using a na ive Bayes model. The accuracy of the classification was 86% using 10-fold cross validation.

In the current paper, to improve the quality assessment for perinatal monitoring, we propose a simpler template-based method using only the DUS signal recorded by a low cost device, thus facilitating its implementation in LMICs.

# 2. METHODS

#### 2.1. Database

The DUS database used in this paper was collected at the John Radcliffe Hospital in Oxford as a part of the study presented in Stroux and Clifford (2014) and Stroux (2016). The study was approved by the NHS Health Research Authority, REC reference: 12/SC/0147 and written consent was obtained from each study subject prior to data collection. Each subject received detailed information on the study protocol and their right to withdraw from the study at any stage of the recording session. This database contained 1D-DUS signals recorded from 17 subjects at a sampling frequency of 44.1 KHz using a hand held transducer (AngelSounds Fetal Heart Detector, Jumper Medical Co., Ltd.) with an ultrasound frequency of 3.3 MHz. Subjects were women with singleton pregnancy, over the age of 18, who were scheduled for a routine CTG. The duration of recording per subject is shown in **Figure 1**.

### 2.2. Segment selection

Each of the 1 min-length DUS signals were labeled by three different annotators who had relevant experience in analyzing cardiac audio recordings for quality. Each reviewer independently labeled each second of the record as good or poor quality. After labeling, each record was split into segments of 3.75 with a 3 s sliding window (i.e., a 0.75 s overlap). The duration of 3.75 was fixed since it is the usual length for computerized analysis of fetal non-stress tests based on the Dawes/Redman criteria (Dawes et al., 1981; Pardey et al., 2002). To ensure that

each segment belongs to only one class, only the segments with all theirsamples of the same class were kept. These segments were assigned to 4 different classes as follows:


A total of 1,990 segments (430 good, 1,062 poor, 292 mostly clean, and 206 mostly noisy quality) were identified. **Figure 2** illustrates the balance of segments across patients. Note that the quality of the recorded signals varies from one patient to another and may change over a single recording session because we observed that for some recordings there are both good and bad quality segments.

The classifier in this work was only trained using poor and good quality segments. The rationale behind this stems from the fact that segments on which one or more experts cannot agree are not meaningful in reporting statistics, since we cannot categorize them in a single class. However, after training the classifiers, the optimal classifier was also evaluated with the mostly clean and mostly noisy segments to determine its capacity for detect intermediate quality segments.

#### 2.3. Preprocessing

The DUS signals were resampled at 4,000 Hz using a leastsquare linear phase anti-aliasing filter. This downsampling does not affect the information content of the signal, since the fetal heart activity corresponds to the DUS signal frequencies below 1,650 Hz for a transducer of 3.3 MHz (Shakespeare et al., 2001). Hence, the Nyquist-Shannon sampling criterion was satisfied after downsampling.

# 2.4. Template-based quality assessment of 1-D Doppler Ultrasound

To assess the quality of the DUS segment, a template-based algorithm was developed. This method consists of 4 stages (**Figure 3**). First, the beats of the DUS signals were estimated using Empirical Mode Decomposition (EMD). Then, using the estimated beats, templates for windows of 15 s were derived. These templates were then optimized in stage 3, and finally, the quality index of the DUS segment was calculated in stage 4. These stages are illustrated in **Figure 3** and explained in the following sections.

This method, as all the remaining methods of this work, were implemented in Matlab and executed in a machine with a standard processor (Intel(R) Xeon(R) CPU E5-2660 v2 @ 2.20GHz).

#### 2.4.1. Beat Detection

Individual cardiac cycles (or beats) were detected using EMD, based on a method presented in Marzbanrad et al. (2014). EMD is an empirical method for decomposing non-stationary and nonlinear signals into a set of components called Intrinsic Mode Functions (IMFs). It is a data-driven method that is able to adapt to the signal properties without requiring a basis function, unlike other time-frequency decomposition methods (Huang, 2014). This characteristic allows EMD to properly analyze non-linear and non-stationary natural processes.

Each extracted IMF satisfies 2 properties: firstly, the number of maxima and minima and the number of zero crossing should differ at most by 1; secondly, the mean value between the envelope of the local maxima and the envelope of the local minima must be zero at any point. To obtain the IMFs, EMD uses an intuitive algorithm called "sifting procedure." It is an iterative procedure, which finds all the IMFs of the signal until the difference between output and the input of the sifting procedure becomes a monotonic function. More details of the method can be found in Huang et al. (1998).

To find the beats from the DUS signals an algorithm was developed to allow switching between the first four IMFs, which were obtained over 4 s windows. For each of these IMFs, the peaks were detected based on the positive first derivative and negative second derivative criteria. Then, using the identified peaks, the IMFs were enveloped to obtain four IMF envelopes for each window. To detect the best envelope for segmentation, a metric based on the standard deviation of the peak to peak intervals (PPIs) was applied. Namely, the IMF with the minimum average of the standard deviation of PPIs was selected as the optimum IMF. To deal with the possible mismatching of the selected IMFs in adjacent windows, a short overlapping window of 1 s was used to correct missing or double identified peaks. The peaks of the optimum IMFs were selected as possible beat locations. These peak locations were further corrected through a moving windows of 5 PPIs, replacing the middle PPI by the average of the rest in the window, if they differed by more than 20% (Marzbanrad et al., 2016a). The corrected peaks were set as beat location, and were used to segment the 1D-DUS signal into Beat to Beat Intervals (BTBI).

Continuous Wavelet transform (CWT) was then applied to the DUS signal over 25 s windows. In this work, the CWT was applied using second order complex Gaussian function as Wavelet mother. Moreover, only the signal decomposed at scale 3 was selected since it was found that the 3rd level is the most relevant for detecting value movement. This scale contained frequencies below 1,000 Hz, which mainly reflect the fetal heart activity; valve motion is around 990 Hz for a transducer of 3.3 MHz (Murata and Martin Jr, 1977), and wall velocities between 257 and 429 Hz for a transducer of 3.3 MHz (Shakespeare et al., 2001).

After applying the CWT, the envelope of the absolute value of the decomposed signal was estimated by interpolating the maxima. This envelope was then smoothed using a low pass filter and segmented into cardiac cycles using the estimated beat locations. Each segment was normalized by subtracting its mean and dividing by its standard deviation. These normalized segment were used to generate the templates for the signal quality assessment.

#### 2.4.2. Initial Template Generation

Using the normalized cardiac cycle segments, the initial templates were calculated using a window of 15 s. The length (L) of the template was calculated as the average of BTBIs in each 15 s window. The initial template of the window was determined by averaging the segments starting at the beat location and lasting at length L. This procedure was repeated for each window, thereby obtaining an initial template for each window of 15 s.

#### 2.4.3. Updating Templates

The initial template of each window was updated based on the correlation function. For each window, the correlation of the template and the segments starting in a beat location and lasting at length L was calculated. The window template was updated by averaging only the segments with a correlation (r) greater than or equal to 0.6. In case the initial template of a window did not have a correlation of r ≥ 0.6 for at least 20% of the beats, the template was assumed as invalid, and replaced by the one from the previous window. If the initial template of the previous window was invalid, the one from the next window was selected.

#### 2.4.4. Signal Quality Metrics

After updating templates for each window, the quality indices were calculated as the correlation of the segments with the template in their corresponding window. The correlations was calculated in four ways:


For all methods, any negative values of these SQI (negative correlation) were set to zero.

# 2.5. Sample Entropy and Power Spectrum Density (PSD)

In addition to the Template-based SQIs, two other key features were estimated from the DUS signals. The first one was sample entropy (Hs), which has shown a promising potential for discriminating between good and poor quality DUS segments (Stroux and Clifford, 2014). Sample entropy measures the regularity of a signal by finding reoccurring patterns in it. To this end, three parameters are defined: the length of the signal N, the pattern length m and the matching tolerance r. Sample entropy is defined as the negative logarithm of the probability that a time series of length N with reoccurring pattern of length m within a set tolerance of r, also has reoccurring patterns of length m + 1 under the same tolerance constraint. In this work, the sample entropy was calculated setting the parameters m = 2, and r as 0.1 times the standard deviation of the input time series. The entropy was calculated using the procedure described in Richman and Moorman (2000).

The second additional feature extracted was the Power Spectrum Density (PSD) ratio. This feature was used in order to the evaluate the power of the DUS signals at different frequency ranges. The range for calculating the ratio was determined using a grid search. Since cardiac movements are associated with a Doppler frequency range of 100 to 600 Hz using a 2 MHz transducer (Wheeler et al., 1987), which translate to a scaled range of 165 to 990 Hz for 3.3 MHz transducer, the ranges of values of the grid search were fixed from 80 to 400 Hz and from 580 to 900 Hz for the low and high frequency interval limits, respectively. For each possible pair of values, the capacity for discriminating between good and poor quality segments was measured using the earth mover's distance. The range with the highest earth mover's distance between the distribution of the ratios of good and poor classes was found to fall in the range 160– 660 Hz. Thus, the PSD ratio of each DUS segment was calculated by dividing the power spectrum contained in the interval [160 − 660 Hz] by the total power, thereby measuring the percentage of power associated with cardiac movements.

#### 2.6. Feature Vectors

Applying the template-based method resulted in four different SQIs for each estimated beat of the segment, thus obtaining a total of 4N<sup>b</sup> indices by segment, where N<sup>b</sup> is the number of beats of the segment. As the number of beats varied for each segment, we selected the median value of each quality index of the segment as the final SQI. Thus, each segment had only one value for SQI1, SQI2, SQI3, and SQI4. Finally, the sample entropy and the PSD ratio were added to the feature vector, thereby obtaining a total of six features for each segment.

#### 2.7. Classification

The above features were used to train an SVM classifier. SVM is a classifier that finds the best hyperplane that maximizes the margin between two classes. When the data are not linearly separable, a kernel function is used to transform the data to a different space in which the data can be separated. In this work, a Gaussian radial basis function kernel was used:

$$k(\mathbf{x}\_i, \mathbf{x}\_j) = \exp\left(\frac{-|\mathbf{x}\_i - \mathbf{x}\_j|^2}{2\sigma^2}\right),\tag{1}$$

where x<sup>i</sup> and x<sup>j</sup> are feature vectors, and σ is a free parameter of the kernel. A large value of σ increases the bias but reduces the variance of the classifier and a small value causes the opposite effect. To find the best value for a given training set, σ is usually tunned using heuristic methods or a brute force search.

The class prediction, y, of a given feature vector, x, is calculated using the dual representation of SVM as:

$$\chi = \text{sgn}\left(\sum\_{i=1}^{N} \alpha\_i \wp\_i k(\chi\_i, \mathfrak{x}) + b\right),\tag{2}$$

where x<sup>i</sup> is the i-th feature vector of the training set and y<sup>i</sup> = [−1, 1] is its class; α ≥ 0 are Lagrange coefficients obtained by quadratic optimization; b is the intercept of the margin; and k(x<sup>i</sup> , xj) is the kernel function (Equation 1). The α coefficient is only greater than 0 for those points that are in the margin. These points are called support vectors. In addition to the parameters of Equation 2, SVM has a hyperparameter called the soft margin constraint (C). This parameter regularizes the margin allowing the cost function to ignore some points to establish an adequate margin for the training set. More details concerning the SVM can be found in Abe (2005). In this work, the SVM parameters C and σ were optimized using five-fold cross-validation and a grid search on the training set. The grid search was defined by C ∈ {2 −3 , 2−<sup>1</sup> , ..., 2<sup>5</sup> } and σ ∈ {2 −5 , 2−<sup>2</sup> , ..., 2<sup>2</sup> }.

#### 2.8. Method performance assessment

To assess the performance of the method proposed here, the dataset was split into two equal subsets; the training and test sets. The the training set was used to determine the combination of features most relevant for assessing the quality of DUS segments. The test set provides an assessment of the accuracy on an independent dataset.

To split the dataset into two equal subsets (training and test sets), the subjects were ranked based on their number of good segments in descending order. Then, the data of each of the subjects were alternately assigned to the subsets. In other words, the first subject's samples were assigned to the training subset, the second ones to the test subset, the third one to the training subset, and so on. As the number of patients was odd (17), the samples of the last subject were assigned to the subset with the lowest number of segments.

The best combination of features was found by calculating the accuracy of all possible feature combinations on the training set. Since the dataset presented an imbalance among classes (**Figure 2**), the accuracy was calculated using stratified five-fold cross validation with bootstrapping. Specifically, the accuracy of each feature combination was determined as follows: subjects of the training set were split into 5 folds. For each fold, 120 signals (60 per class) samples from the subjects of the fold were randomly selected using sampling with replacement (bootstrapping). The selection was performed in proportion to the subjects' sample quantity in each fold. The rationale behind this validation process is that the bootstrap applied to the cross validation folds adjusts the class imbalances, which is a critical factor for SVM classifiers (Chawla et al., 2004). Moreover, as the cross validation did not assign samples of the same subject to different folds, it provided an unbiased accuracy estimation. To obtain a more reliable accuracy, the described validation process was repeated 100 times, assigning subjects into different folds at each repetition.

For each iteration of the five-fold cross validation, the training set was normalized by subtracting the mean of the respective feature vector and dividing by its standard deviation. The test set was normalized using both mean and standard deviation derived from the training data. The cross validation accuracy of each iteration was averaged by selecting the median of the 5 accuracy values of the folds. Likewise, the accuracy of the 100 repetitions was selected as the median of the 100 accuracy values. This procedure was performed for each combination of features.

In addition to the overall accuracy of the classifiers, the sensitivity and specificity were also estimated. Sensitivity was defined as the proportion of good quality segments properly classified, whereas specificity denoted the portion of poor quality segments correctly classified.

To determine the capacity of the method for predicting intermediate quality segments, a SVM classifier was trained with the good and poor segments of the test set using the most common parameters C and σ for the 100 bootstrap repetitions. Once the best combination of features was determined (maximizing accuracy, then specificity), the classifier was fixed and assessed on the test using the same bootstrap cross-validation validation procedure used for the training set. Finally, the probability of belonging to good class was also estimated for the mostly noisy, and mostly clean segments without retraining to assess the performance of the classifier on all data.

#### 3. RESULTS

#### 3.1. Feature Selection

**Table 1** presents the median accuracy, sensitivity and specificity of the best combination of input features for up to 6 possible features. As can be seen, the classifier was able to classify the quality of a DUS segments with up to 85.8% accuracy using either the combination SQI2 and sample Entropy (HS), or the combination SQI2, PSD ratio and HS. The accuracy of these two combinations of features resulted in a statistically significant improvement over the use of only one feature, H<sup>S</sup> (p <0.05, one-sided Wilcoxon rank-sum test). The results for all possible combination of features are presented in the appendix (**Table A1**).

It can also be seen from **Table 1** that the sensitivity tended to decrease with an increase in the number of features, whereas specificity steadily increase until three features were used. Since the combination of both two and three features leads to equivalent accuracy, the combination SQI2, PSD ratio and H<sup>S</sup> was chosen to be evaluated on the test set, since this maximizes specificity, and reduces the chances that a poor quality segment is labeled as good quality.

#### 3.2. Test Set Performance

**Table 2** displays the accuracy, sensitivity and specificity of the combination of the SQI2, PSD ratio and H<sup>S</sup> features. The median accuracy of this classifier using this combination was similar to the highest median accuracy achieved on the training set. However, the interquartile range for the test set was almost twice than that for the training set, indicating that the test set may exhibit a higher heterogeneity of features. Both sensitivity and specificity exceeded 90%.

### 3.3. Performance of classifier on intermediate quality segments

The classes of the mostly clean and mostly noise segments of the test set were also predicted using the same classifier [using



IQR indicates inter-quartile range; ‡ indicates a significant improvement (Wilcoxon ranksum test, P<0.05) of a given feature combination compared to using a combination with one less feature.

TABLE 2 | Performance of the classifier averaged over 100 five-fold cross validation runs balanced with bootstrapping for the test set (with ) using SQI2, PSD ratio, and sample entropy (SQI2,PSD,Hs) as features.


SQI2, PSD ratio and sample entropy (HS) as features]. Note that these segments were not used in training. **Figure 4** shows the relative distribution of output probabilities from the classifier of belonging to the good quality class for all four types of segments. The classifier established a probability threshold of 0.5575 for distinguishing between good quality and poor quality segments. The percentage of segments which lay above the threshold for good quality and mostly clean segments were 86.53 and 69.06%, respectively. On the other hand, the percentage of segments that lay below the threshold for poor quality and mostly noise segments were 96.50 and 63.69%, respectively.

## 4. DISCUSSION

The results presented here suggest that it is possible to accurately classify the DUS quality using SQIs derived from DUS signals alone. Among the extracted features, sample entropy and PSD ratio provided suitable discrimination between good and poor quality segments, which is consistent with previous works (Stroux and Clifford, 2014; Marzbanrad et al., 2015). However, the addition of our proposed template-based method, particularly after linear resampling of the beats to match the running template (SQI2), provided a statistically significant improvement in accuracy (see **Table 1**). Either combinations SQI2 and H<sup>s</sup> or SQI2, PSD ratio and H<sup>s</sup> resulted in a statistically significant accuracy; nevertheless, in order to maximize specificity, the combination of SQI2, PSD ratio and H<sup>s</sup> was selected for assessing the classifier on the test set.

The selected features achieved an accuracy of 85.8% on the test set, thus suggesting that this metric is suitable for quality assessment based only on DUS signals. Although this feature combination exhibited more variance on the test set than on the train set, the achieved accuracy indicates that the model was not overfitted, and its complexity of three features is viable for assessing DUS quality. Furthermore, the balance toward specificity provided by the three chosen features (SQI2, PSD

FIGURE 4 | Distributions of classifier probability outputs for DUS segments of test set for each of the four classes (n.u. stands for normalized units). The threshold of belonging to the Good class was fixed at 0.56 for the classifier. The majority of the distribution of the Good and Mostly Clean classes lies above this threshold, whereas the majority of Poor and Mostly Noise classes lies below this threshold, as was expected. The probability distributions were smoothed using a normal kernel function (Bowman and Azzalini, 1997).

ratio, and sample entropy) ensures a high number of good quality segments is preserved, as well as small number of false-positive segments.

The best combination of features also showed an adequate capacity for classifying segments associated with intermediate quality zones (mostly noisy and moistly clean segments). Although both mostly clean segments and mostly noisy segments exhibited a mostly flat distribution, their centers where more closer to the good quality center and poor quality centers, respectively as it was expected. Specifically, almost 70% of the mostly clean probabilities laid above the SVM prediction threshold, whereas higher than 63% of the mostly noise probabilities laid below the SVM prediction threshold. This discrimination ability for the two ambiguous classes indicates the potential of the approach outlined in this work.

Regarding the template-based SQIs, the EMD-based approach appeared to facilitate identification of beat intervals, since the correlation with each template was generally high. Small offsets in the relative start and end point of each beat were mitigated by the use of resampling prior to correlation. The segmentation facilitated beat-by-beat quality assessment, which is the first step toward detecting fetal abnormalities from DUS signals. The template optimization process obtained representative templates for quality assessment since the initial template was only averaged with those segments which exhibited a moderate or strong correlation (r ≥ 0.6) with the initial template (average of first N beats).

Although Stroux and Clifford reported a higher accuracy (95.14%) on the same database (Stroux and Clifford, 2014), their work cannot directly compared to the current work since different statistical validation approaches were used. Specifically, Stroux and Clifford trained on two thirds of the data set and held out one third for testing, with no cross validation. In this work, stratified five-fold cross validation was used with bootstrapping (repeated 100 times), with subject stratification across different folds in each repetition. The accuracy obtained in this work cannot directly compared to that of Marzbanrad et al. (2015) since they tested their method with a different dataset. However, our method can be compared to the aforementioned previous works by analyzing the effect of adding the index SQI2 to common features of the other works, namely, sample entropy. As was previously showed in **Table 1**, by using the SQI2 feature in addition to sample entropy, the accuracy statically significantly increased.

Another advantage of our method over previous works is that the proposed method does not need additional sources, such as accelerometer data (Stroux and Clifford, 2014) or an fECG signal (Marzbanrad et al., 2015) to assess the DUS quality. A key advantage of using only DUS signal is that the recording process is simple, facilitating the use of this technology by non-experts in low-resource settings. Finally, using only one source for quality assessment reduces health screening costs, facilitating its use in LMICs.

Despite the promising results, one limitation of the current method is that it was only tested using DUS signals recorded by professional midwives in hospital settings. The LMIC context is often severely resource constrained and there is a lack of widespread training for midwives, particularly in the use of technology. Consequently, signal quality is likely to be lower in recordings taken in LMICs. The noise content may also be different if the audio cable is not incorrectly inserted, introducing ambient sounds such as animals, extreme weather, and interference from non-hospital electronics. Nevertheless, the template-based method proposed here could be adjusted to specific conditions with a relevant training set.

Another possible limitation may be that the introduced method was only tested using one database labeled by three annotators. As DUS quality annotation is prone to interobserver variability, testing the method with datasets annotated by different experts may reduce the accuracy. However, the high accuracy achieved by the combination of sample entropy, PSD ratio and SQI2 used in this work, provides optimism for the use of the template-based method for different datasets, especially with retraining.

# 5. CONCLUSION

The work presented in this article proposed a template-based method to assess the quality of 1D-DUS signals recorded by a low cost device. The introduced template-based indices provided a simpler method based on only DUS signals, thus facilitating its implementation in LMICs. The approach described in this work can provide the operator with an accurate and timely feedback on the quality of the recordings, to allow discarding the low quality signals in real time and prompt users to rerecord data. Therefore, this quality assessment technique could potentially facilitate reliable fetal monitoring by non-experts toward reducing perinatal health burdens.

# AUTHOR CONTRIBUTIONS

CV implemented the second version of the method, performed experiments, analyzed results, wrote the paper. FM designed the method, implemented the first version of the method, she reviewed edited the paper. LS recorded the One Doppler Ultrasound signals used in this work. GC is the supervisor of CV; he guided the experiments performed by CV; he reviewed and edited the paper.

# ACKNOWLEDGMENTS

CV acknowledges support through a Fulbright Scholarship. LS acknowledges the support of the RCUK Digital Economy Programme grant number EP/G036861/1 (Oxford Centre for Doctoral Training in Healthcare Innovation) and of the Oxford Centre for Affordable Healthcare Technology. GC acknowledges the support of the National Institutes of Health, the Fogarty International Center and the Eunice Kennedy Shriver National Institute of Child Health and Human Development, grant number 1R21HD084114-01 ("Mobile Health Intervention to Improve Perinatal Continuum of Care in Guatemala"). Sponsorship was also received from ARM through the Center for Affordable Healthcare Technology at Kellogg College, Oxford for early data collection and prototyping.

#### REFERENCES


from doppler ultrasound signals without using fetal electrocardiography," in Computing in Cardiology Conference (CinC), 2014, (Cambridge, MA: IEEE), 485–488.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Valderrama, Marzbanrad, Stroux and Clifford. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# APPENDIX

TABLE A1 | Median classification performance of the 100 five-fold cross validation balanced with bootstrapping for all the possible feature combinations.


#### TABLE A1 | Continued


The table is grouped for feature vectors of the same length. For each combination of features, the median and interquartile range of the accuracy rate of the 100 repetitions are shown.

(Continued)

# **Fetal Heart Sounds Detection Using Wavelet Transform and Fractal Dimension**

*Elisavet Koutsiana<sup>1</sup> \*, Leontios J. Hadjileontiadis 2,3 \*, Ioanna Chouvarda<sup>1</sup> and Ahsan H. Khandoker 3,4*

*1 Laboratory of Medical Informatics, The Medical School, Aristotle University of Thessaloniki, Thessaloniki, Greece, <sup>2</sup>Department of Electrical and Computer Engineering, Aristotle University of Thessaloniki, Thessaloniki, Greece, <sup>3</sup>Department of Electrical and Computer Engineering, Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates, <sup>4</sup>Department of Biomedical Engineering, Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates*

Phonocardiography is a non-invasive technique for the detection of fetal heart sounds (fHSs). In this study, analysis of fetal phonocardiograph (fPCG) signals, in order to achieve fetal heartbeat segmentation, is proposed. The proposed approach (namely WT–FD) is a wavelet transform (WT)-based method that combines fractal dimension (FD) analysis in the WT domain for the extraction of fHSs from the underlying noise. Its adoption in this field stems from its successful use in the fields of lung and bowel sounds de-noising analysis. The efficiency of the WT–FD method in fHS extraction has been evaluated with 19 simulated fHS signals, created for the present study, with additive noise up to (3 dB), along with the simulated fPCGs database available at PhysioBank. Results have shown promising performance in the identification of the correct location and morphology of the fHSs, reaching an overall accuracy of 89% justifying the efficacy of the method. The WT–FD approach effectively extracts the fHS signals from the noisy background, paving the way for testing it in real fHSs and clearly contributing to better evaluation of the fetal heart functionality.

**Keywords: fetal heart rate, fetal heart sound, fetal phonocardiogram, wavelet transform, fractal dimension thresholding**

# **INTRODUCTION**

Fetal heart rate (fHR) observation is important for proper fetal well-being assessment during the period of pregnancy. Electronic fetal monitoring (EFM) is a significant tool for the obstetrician, in order to perform various tests at different stages of gestation to estimate the fetal health. The typical examination until the 28th week of pregnancy is composed of continuous measurements of the fetus growth, while at the stage of 29–40 weeks, the monitoring of fetal movement, fetal respiration, fHR, and others (Adithya et al., 2017) are included. In current practice, the examination of the fetus is performed by means of ultrasonic-based equipment such as Doppler ultrasound and cardiotocogram (Nassit and Berbia, 2015).

Although Doppler ultrasound and cardiotocogram are the typical fetal observation devices, these techniques present some limitations, mainly because of the cost of the monitoring devices and the complexity of their use, demanding an expert during data acquisition. Moreover, it has not been established that the frequent and long-term exposure to ultrasound energy has no effect on either the fetus or the mother (Salvesen, 2002).

Existing standards of fetal monitoring estimate the fetus and the mother physiology with repetitive examinations. However, complications may occur during pregnancy, i.e., fetal deaths,

#### *Edited by:*

*Rajat Mittal, Johns Hopkins University, United States*

#### *Reviewed by:*

*Satoru Ikenoue, Keio University, Japan William Reid Thompson, Johns Hopkins School of Medicine, United States*

#### *\*Correspondence:*

*Elisavet Koutsiana koutsiana@auth.gr; Leontios J. Hadjileontiadis leontios@auth.gr*

#### *Specialty section:*

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Bioengineering and Biotechnology*

> *Received: 20 March 2017 Accepted: 03 August 2017 Published: 08 September 2017*

#### *Citation:*

*Koutsiana E, Hadjileontiadis LJ, Chouvarda I and Khandoker AH (2017) Fetal Heart Sounds Detection Using Wavelet Transform and Fractal Dimension. Front. Bioeng. Biotechnol. 5:49. doi: 10.3389/fbioe.2017.00049*

**27**

preterm delivery, hypoxia, and other, which have no specific prevision. Even though the literature is not robust about the risks and their relation to the EFM, long-term fHR monitoring has proven to be an effective approach for better accuracy in the clinical examination of the fetus (Martin, 1998).

A passive alternative for long-term monitoring of the fetus is the fetal heart auscultation. It is a non-invasive method that records the vibroacoustic signals from the abdominal surface. The acoustic signal produced by the fetal heart sound (fHS) can be visually depicted in the fetal phonocardiograph (fPCG). We can separate the fetal heartbeat into two sub-beats, the systolic beat S1 and the diastolic beat S2, which follows S1. The S1 and the S2 subbeats are generated by the vibratory components of the fetal heart valves closure. The S2 sub-beats present smoother morphology than the S1, making harder the detection of their location. A heart cycle consists of the S1 and S2 sub-beats.

The research of fPCG signals aims to segment the S1 and S2 sub-beats, in order to study the wavelet morphology of the fHSs and the fHR variability. Long-term monitoring of the fHSs reveals information about the fetus growth and functionality. Although there is not enough knowledge about fHS morphology, in order to indicate any pathological conditions, the study of the fPCG signals have shown promising results to the extension of the EFM and the physical examination of the fetus (Adithya et al., 2017).

Auscultation is a low-cost and non-invasive method as it captures the acoustic signal of the fHSs. Moreover, the phonocardiogram device is a flexible method that does not need an expert to record the signals. The mother can take long-term recordings during the day or night and afterward, the doctor can examine the signals and have a more complete overview of the fetus functionality.

Nevertheless, fetal auscultation has many challenges. Because of the place of the fetus in the maternal abdominal, the fPCG signals are loaded with noise from various sources such as maternal heart sounds, digestive sounds, maternal and fetus respiration movements, external noise, and others (Várady et al., 2003; Cesarelli et al., 2012). In the noisy fPCG signals, the fetal heartbeats are often masked by other components, consequently it is difficult to detect without applying robust signal processing methods.

Throughout the years, various signal processing approaches for de-noising the fPCG signal have been examined and proposed (Unser and Aldroubi, 1996; Messer et al., 2001; Várady et al., 2003; Xiu-Min and Gui-Tao, 2009; Chourasia and Mittra, 2010; Chourasia et al., 2011, 2014). Among them, Khadra et al. (1991) were the first to suggest the wavelet transform (WT) as a useful tool for the analysis of heart sounds. Following, many researchers concentrated on the study of wavelet-based techniques for these signals. Vaisman et al. (2012) proposed the WT as a de-noising tool for the determination of the fHR. At the same time, Kovács et al. (2011) used autocorrelation technique, WT, and matching pursuit for the evaluation of fHS. Recently, Chourasia and Tiwari (2013) designed a new wavelet basis function for de-noising the fPCG signals.

The present study was motivated from a previously proposed method of Hadjileontiadis for the separation of lung and bowel sounds from the background noise (Hadjileontiadis, 2005). The latter technique uses a scheme of WT for de-noising the signals and also fractal dimension (FD) analysis for the detection of lung and bowel sounds. The so-called WT–FD filter introduces an alternative way to the enhancement of bioacoustic signals, applicable to any separation problem involving non-stationary transient signals mixed with uncorrelated stationary background noise (Hadjileontiadis, 2005).

In this study, the WT–FD method is suggested for the case of fPCG signals, to effectively locate and extract the fetal heartbeat from the underlying noise. Due to the highly noisy environment and the low acoustic energy of the fetal heartbeat, WT is an efficient method that decomposes the signal into multiple levels for the subtraction of the unwanted stationary noise. Moreover, the method is flexible since it uses short windows at high frequencies and long windows at low frequencies making the wavelet function more similar to the waveforms of the signal. Furthermore, FD analysis is frequently used in biomedical signal processing. There are studies of FD performance at electroencephalograms for the detection of the onset of epileptic seizures and also at electrocardiogram signals for the classification of arrhythmia with satisfying results (Mishra and Raghav, 2010; Polychronaki et al., 2010).

The rest of the paper is formed as follows. Section "Mathematical Background" describes the mathematical background of WT and FD definitions, while Section "The WT-FD Method" presents the proposed method. Section "Implementation and Evaluation Issues" describes the databases that the method was tested and the general indices that used for its evaluation. Finally, Section "Results" confers some experimental results, which evaluate the efficiency of WT-FD algorithm in fPCG signals, and Section "Concluding Remarks" concludes the paper with suggestions for future work.

# **MATHEMATICAL BACKGROUND**

#### **Wavelet Transform**

Wavelets are families of functions ψ*a*,*b*(*t*) generated from a singlebase wavelet ψ(*t*) called the "mother wavelet," by dilations and translations (Hadjileontiadis and Panas, 1998; Olkkonen, 2011), i.e.,

$$
\Psi\_{a,b}(t) = \frac{1}{\sqrt{a}} \Psi\left(\frac{t-b}{a}\right), \ a > 0, \ b \in \mathbb{R}, \tag{1}
$$

where *a* is the dilation (scale) parameter and *b* is the translation parameter.

In the past few decades, wavelet analysis has been proved to be an important tool in biomedical engineering. The use of WT in fPCG signals is driven by the nature of the signals itself. Explosive peaks in the time domain produce large coefficients over the wavelet scales, while the noisy background dies out swiftly with increasing scale. In WT, the signal is decomposed into coarse and detail information using a pair of finite impulse response filters (and their adjoins), which are low-pass and high-pass, respectively (Hadjileontiadis, 2005). The process can be described as a tree, which at each step decomposes the low-pass filter into further lower and higher frequency coefficients. Thus, the original signal is decomposed into coefficients of lower resolution, and the high frequency coefficients are not analyzed any further. This scheme is a wavelet-based multiresolution decomposition, and it is known Koutsiana et al. fHSs Detection

as Mallat algorithm (Mallat and Peyré, 2009). The procedure that uses the coarse and the detail coefficients and yields back to the original signal is multiresolution reconstruction.

In the proposed de-noising method, the decomposition– reconstruction scheme was based on the orthonormal bases and the quadrature mirror filters introduced by Daubechies (1988). This wavelet family was chosen because of the morphology of the mother wavelet, comparatively the waveforms of the fPCG and the testing of other wavelet families.

#### **Fractal Dimension**

Fractals are mathematical sets, which describe many natural phenomena with geometrical complexity (Mandelbrot, 1982; Esteller et al., 2001). The term "fractal dimension" can more generally refer to any of the dimensions commonly used for fractals characterization (e.g., capacity dimension, correlation dimension, information dimension, Lyapunov dimension, and Minkowski–Bouligand dimension) (Hadjileontiadis, 2005). More accurately, the FD is a priceless tool that reflects the signal complexity in the time domain. Here, FD was adopted as a means to detect the most important WT coefficients that correspond to the fetal heartbeat in the WT domain, resulting, simultaneously, in significant computational savings.

The FD technique is performed using a sliding window of *W* = int(0.05*·Fs*) samples length, where int(*·*) indicates the integer part of the argument, the constant is empirically set at 0.05 justifying the efficient performance of the algorithm, and *F<sup>s</sup>* denotes the sampling frequency of the signal. It is noticed that when the *W* window is small, too many false FD peaks are generated and when it is big, the estimated FD is smoothed so the algorithm chooses the false peaks.

Let the processing signal be an *N*-sample vector. Then, the *W*sample window is one-sample shifted along the *N*-sample input vector in order to obtain point-to-point values of the estimated FD. Every estimated FD obtained with the sliding window is assigned to its midpoint. In this way, the length of the final sequence of the FD(*i*) is lower than *N*. This length is extended to comply with the *N*-sample length of the original input vector, assigning the FD(1) and FD(*N − W* + 1) estimated values to the first and last half of the *W −* 1 missing values, respectively. In this study, we used the Katz's definition of FD as it is proposed by Hadjileontiadis and Rekanos (2003) for the detection of explosive lung and bowel sounds.

According to Katz (1988), the FD of a curve defined by a sequence of *N* points is estimated by

$$\text{FD} = \frac{\log\_{10}(n)}{\log\_{10}\left(\frac{d}{L\_{\epsilon}}\right) + \log\_{10}(n)},\tag{2}$$

where *L*<sup>c</sup> is the total length of the curve, realized as the sum of distances between successive points, i.e.,

$$L\_{\mathfrak{c}} = \sum\_{i=1}^{N-1} \text{dist}(i, i+1),\tag{3}$$

where dist(*i*,*j*) is the distance between the *i* and *j* points of the curve; *d* is the diameter estimated as

$$d = \max[\text{dist}(i, j)], \qquad i \neq j, \qquad i, j \in [1, N], \tag{4}$$

for curves that do not cross themselves; usually, the *d* diameter is estimated as the distance between the first point of the sequence and the point of the sequence that provides the farthest distance, i.e.,

$$d = \max[\text{dist}(1, i)], \qquad i, j \in [2, N], \tag{5}$$

and *n*<sup>s</sup> is the number of steps in the curve, defined as

$$n\_{\mathfrak{s}} = \frac{L\_{\mathfrak{c}}}{\alpha},\tag{6}$$

where α denotes the average step, i.e., the average distance between successive points.

### **THE WT–FD METHOD**

#### **WT–FD Iterative Procedure**

The WT–FD method is an iterative procedure performed in order to achieve the best separation of fetal heartbeat from the superimposed noise. The amplitude normalized *N*-sample input vector *X*[*n*] (*n* = 1, *. . .*, *N*), is subjected to the WT–FD technique and is separated into two parts, i.e., *X k S* [*n*] and *X k <sup>U</sup>* [*n*], the nonstationary desired signal and the stationary background noise, respectively. After that, the process continues iteratively with the vector *X k <sup>U</sup>* [*n*] serving as a new *X*[*n*] input signal to the next iteration, and the resulted vectors across all *L* iterations, i.e., *X*[*n*] 1:*L* , are used for the final reconstruction (Hadjileontiadis, 2005). The iterative procedure stops when the following stopping criterion is satisfied:

$$\text{STC} = \left| E\left\{ X\_U^{k-1}[n]^2 \right\} - E\left\{ X\_U^k[n]^2 \right\} \right| < \varepsilon,$$

$$k\text{-th iteration, } n = 1, \dots, N,\tag{7}$$

where *E*{*·*} denotes the expected value. The parameter ε is a small positive number (0 *<* ε *≪* 1*.*0) that corresponds to the desired accuracy in procedure. The initial value of *X* 0 *<sup>U</sup>* [*n*] is considered to be equal to 0. When the STC criterion is satisfied after *L* iterations, the final reconstruction of the signal is achieved with the *X k S* [*n*] vectors as follows:

$$X\_{\text{REC}}[n] = \sum\_{k=1}^{L} X\_{\text{S}}^{k}[n], \ n = 1, \ldots, N,\tag{8}$$

A schematic representation and further details about the WT–FD filter can be found in Hadjileontiadis (2005).

#### **WT Coefficient Estimation and Selection**

In this study, the Daubechies 4 wavelet family (Daubechies, 1988) has been chosen for de-noising the signal. As described in Section "Wavelet Transform," WT decomposes the input fPCG signal *X*[*n*] (*n* = 1, *. . .*, *N*) into *R* detail coefficients WT*m*[*n*] (*m* = 1, *. . .*, *R*). The number *R* of the adjustment resolution scales is estimated by log<sup>2</sup> *N*. An example of an fPCG signal is presented in **Figures 1A–H**, where the original signal is decomposed into seven levels. It is clearly depicted that the first WT level contains only noise and the last three do not contain any important components of the signal. Hence, from the *R* estimated coefficients

the algorithm selects those including important information and leaves out those including background noise, as described next.

First, from the *R* estimatedWT resolution levels, the first *D* ones are discarded according to the following criterion:

$$D = \min \{ \lambda : \mathfrak{n}'\_1 - \mathfrak{n}'\_{\lambda+1} \le 0.4 \}, \tag{9}$$

and from the *J* = (*R − D*) coefficients, the first *M* ones are selected according to the following criterion (Hadjileontiadis, 2007):

$$M = \min\{\lambda : |\mathfrak{n}'\_{\lambda}| > p^{\wedge} \, |\mathfrak{n}'\_{\lambda+1}| \le p^{\wedge} \, \mathfrak{n}''\_{\lambda} > 0\},\tag{10}$$

with

$$\mathfrak{m}\_{\lambda} = 1 - \frac{\sum\_{i=1}^{\lambda} E\{\mathrm{WT}\_{i}(n)^{2}\}}{\sum\_{i=1}^{J} E\{\mathrm{WT}\_{i}(n)^{2}\}}, \ \lambda = 1, 2, \dots, J \ \mathfrak{n} = 1, 2, \dots, N,\tag{11}$$

where η *′* <sup>λ</sup> and η *′′* <sup>λ</sup> denote, respectively, the first and second derivatives of η<sup>λ</sup> with respect to λ, *p* is a small number close to 0 that serves as a threshold, which accounts for the fluctuation of the first derivative around 0, and *E*{*·*} denotes the expected value; here, *p* was empirically set equal to 0.01.

#### **FD-Based S1 and S2 Selection**

The fHS segmentation is performing using the FD method across the selected WT*j*[*n*] (*j* = 1, *. . .*, *M*) WT level. Specifically, the windowing Katz definition of FD as it is described in Section "Fractal Dimension" is performed at every selected coefficient. Then, the estimated FD*<sup>i</sup> j* [*n*] (*i*-th iteration*, j*-th selected coefficient) are fed to the FD-peak peeling algorithm (FD-PPA), as it is proposed by Hadjileontiadis (2005), in order to automatically detect the FD*<sup>i</sup> j* peaks. Through a self-adjusted iterative procedure, the FD-PPA iteratively "peels" the estimated FD signal, gradually gathering those parts that construct its peaks, resulting in the FDPP*<sup>i</sup> j* [*n*] sequence as it is shown in **Figures 2A–C**. Hence, the algorithm aims to search for the lower peaks, such as the S2 fetus heartbeats, which correspond to the low amplitude coefficients.

In the present study, each WT*j*[*n*] (*j* = 1, *. . .*, *M*) coefficient is separated in smaller epochs for better FD assessment. Therefore, the FD estimation is more accurate considering the lower peaks. The normal duration of a fetal heart cycle is 430 ms and, consequently, a mean value for each epoch is at 430 ms in order to contain at least one S1 and one S2 heartbeat. Every epoch is fed in the FD-PPA iteration procedure and then reunited in the WT*j*[*n*] (*j* = 1, *. . .*, *M*) coefficient estimation.

The FD-PPA iteration procedure starts with a threshold operation based on the SD of the vector FD*<sup>i</sup> j* [*n*] as follows:

$$\text{pFD}\_{j}^{i} = \begin{cases} \text{FD}\_{j}^{i}, & \text{FD}\_{j}^{i} > \mu^{i} + \sigma^{i} \\ 1.0, & \text{elsewhere} \end{cases}, i = 1, \dots, L\_{1}; j = 1, \dots, M, \tag{12}$$

where μ *<sup>i</sup>* = mean(FD*<sup>i</sup> <sup>j</sup>*) is the mean value of the FD*<sup>i</sup> <sup>j</sup>* vector, σ *<sup>i</sup>* = std(FD*<sup>i</sup> <sup>j</sup>*) is the SD of the FD*<sup>i</sup> <sup>j</sup>* vector, and *L*<sup>1</sup> is the number of the self-adjusted iterations. Thus, the vector *z <sup>i</sup>* = FD*<sup>i</sup> <sup>j</sup> −* pFD*<sup>i</sup> <sup>j</sup>* + μ *i* is created, and 1.0 is the minimum value of the estimated FD sequence. The iterative procedure stops when the following stopping criterion is satisfied:

$$\text{SC}^{i} = \left| E\left\{ \left( z^{i} \right)^{2} \right\} - E\left\{ \left( z^{i-1} \right)^{2} \right\} \right| < \text{acc}, \ i = 1, \ldots, L\_{\text{l}}, \tag{13}$$

where *E*{*·*}, as in the former stopping criterion, denotes the expected value, the parameter acc is a small positive number (0 *<* acc*≪*1.0) that corresponds to the desired accuracy in the procedure, and the initial value of *z* 0 is equal to 0. When the stopping criterion is not satisfied, the vector FD*<sup>i</sup> j* is replaced by the vector *z i* , and it continues the iterative procedure. When the stopping criterion is satisfied, the FD-PPA generates the FDPP*<sup>k</sup> j* sequence of the *j*-th WT coefficient as follows:

$$\text{FDPP}\_{j}^{k} = \sum\_{i=1}^{L\_{1}} \text{pFD}^{i} - (L\_{1} - 1), \; k = 1, \dots, L; \; j = 1, \dots, M,\tag{14}$$

where *L* is the iteration number of the procedure that is described in Section "WT–FD Iterative Procedure." After the FD-PPA implementation, the small peaks that do not correspond to any sound and their duration is less than int(0.015*Fs*), and their normalized amplitude less than 0.25 are removed. Again, int(*·*) indicates the integer part of the argument, the constant is empirically set at 0.015, and *F<sup>s</sup>* denotes the sampling frequency of the signal. Subsequently, the FDPP*<sup>k</sup> j* sequence is generated and thereafter two binary thresholds are constructed, as shown in **Figures 2D,E**. The first binary threshold, i.e., SBTH*<sup>k</sup> j* is used for segmenting the WT coefficients that are related to the desired signal, while the second one, i.e., NBTH*<sup>k</sup> j* is used for segmenting the WT coefficients that are related to the background noise. These two binary thresholds are defined as follows:

3

, the noise binary threshold.

the signal binary threshold, and **(E)** NBTH<sup>1</sup>

$$\text{SBTH}\_{j}^{k} = \begin{cases} 1, \text{FDPP}\_{j}^{k} \neq 1\\ 0, \text{FDPP}\_{j}^{k} = 1 \end{cases},\tag{15}$$

$$\text{NBTH}\_{j}^{k} = \left[1 - \text{SBTH}\_{j}^{k}\right], \ k = 1, \ldots, L; \ j = 1, \ldots, M,\tag{16}$$

The multiplication of the SBTH*<sup>k</sup> <sup>j</sup>* with the WT coefficient gives a set of de-noised signals that create the *X k S* [*n*] vectors as defined in Section "WT–FD Iterative Procedure," while the multiplication of the NBTH*<sup>k</sup> <sup>j</sup>* with the WT coefficients gives the set of the *X k <sup>U</sup>* [*n*]. **Figures 2A–E** gives an example where a working scheme of the proposed method is presented on the third WT level of an input signal. It shows that the FD method successfully detects the location of the sounds by using the binary sequences, and it separates the non-stationary bioacoustics signal from the stationary background noise.

In this study, the final goal is to segment the fHS and separate the S1 from the S2 beats. The decision between S1 and S2 is based on the fact that in a cardiac cycle the diastolic duration is longer than the systolic one (Papadaniil and Hadjileontiadis, 2014). For that reason, the following inequality is checked:

$$\mathcal{S}(2i+1) - \mathcal{S}(2i) < \mathcal{S}(2i+3) - \mathcal{S}(2i+2),\tag{17}$$

where *S*(*l*) is a vector that is created by the binary threshold SBTH*<sup>k</sup> j ,* and it contains the locations of the start and the end of every fetal heartbeat. Moreover, *i* = 1, *. . .*, (*N*1/2) *−* 2, where *N*<sup>1</sup> is the length of the *S*(*l*) vector. If Eq. 17 is true, the interval [*S*(2*i* + 1):*S*(2*i* + 2)] corresponds to S2, otherwise, it corresponds to S1. The first and the last heartbeat of the signal are not determined from this inequality. Hence, they need to be separately defending. For *i* = 1, if Eq. 17 is true, then the second sound [*S*(4):*S*(5)] is S2 and the first sound [*S*(1):*S*(2)] is S1. Respectively, for *i* = (*N*1/2) *−* 2, if Eq. 17 is true, the last sound [*S*(*N*<sup>1</sup> *−* 1):*S*(*N*1)] is defined as S1.

A criterion of each estimated fetal heartbeat amplitude and the distance between fetal heart cycles is also considered for better decision between S1 and S2 beat. In the literature, the mean amplitude of a fetal S1 beat is equal to 0.7 (Cesarelli et al., 2012), and the distance between fetal heart cycles, i.e., between S2 and the following S1, depends on the fHR. The smaller distance between fetal heart cycles is in case of tachycardia and is about 140 ms. Thus, for the decision between S1 and S2 beat, the S1 estimated beat must surpass the 0.5 normalized amplitude and the S2, S1 inter-distance must be outdistance within 130 ms.

#### **IMPLEMENTATION AND EVALUATION ISSUES**

The analysis of this study was applied on a personal computer using Matlab R2015a and tested on simulated databases. Every input signal was tested for 10 s considering *F<sup>s</sup>* = 1,000 Hz, i.e., 10,000 samples.

For the purposes of this research and the algorithmic development of the WT–FD method, a database with fPCG signals was created. Each signal contains simulated S1 and S2 auscultation sounds created by Hadjileontiadis using the model of Chen et al. (1997) and Xu et al. (2001) and adjusted to the duration of fetal heartbeat. The inter-distance between S1 and S2 heart sounds is given by the expression SSID = 210 *−* 0.5 *·*fHR according to Kovacs et al. (2000). Moreover, in order to represent the noise presence, additive white Gaussian noise was used, resulting in signal-to-noise-ratio (SNR) within the range of SNR = [8, 3] dB. The SNR values were computed according to the following steps; measure the power of the signal (*Ps*), convert the given SNR in decibels (SNRdB) to linear scale according to SNRlinear = 10 SNRdB <sup>10</sup> *,* and finally create the noise vector from Gaussian distribution of specific noise variance according to noise = √ *<sup>P</sup><sup>S</sup>* SNRlinear *·* random*,* where random is a vector of normally distributed random numbers with the signal length.

The database consists of signals with different heart conditions corresponding to cases such as tachycardia, bradycardia, and arrhythmia. Specifically, after the 20th week of gestation, the fHR is stabilized between the 110 and the 160 bpm. Thus, for the normal heartbeat signals, the fHR was set at 140 bpm, for the bradycardia signals at 110 bpm, for the tachycardia signals at 180 bpm, and for cases of arrhythmia a range of 80–200 bpm was considered. Hence, many signals with different conditions and different values of additive white Gaussian noise were created and used for testing the present study.

Furthermore, for better assessment of the WT–FD technique, the method was tested on the simulated fPCGs database available *via* PhysioBank.<sup>1</sup> PhysioBank is a large archive of digital recordings of physiological signals and related data for use by the biomedical research community (Goldberger et al., 2000). The simulated fPCG database was created by Cesarelli et al. (2012) and Ruffo et al. (2010). This data set is a series of synthetic fPCG signals related to different fetal states and recording conditions. Simulated fPCG were generated as a sequence of frames, each of which includes simulated S1 and S2 signals, corrupted by noise. These signals are qualified by a range of SNR values that were computed in decibels according to the following formula:

$$\text{SNR} = 10 \log\_{10} \left( \frac{P\_s}{P\_n} \right), \tag{18}$$

where *P<sup>s</sup>* and *P<sup>n</sup>* are the power of fHS and the power of the noise, respectively. The noise source was simulated by generating maternal and fetal noise, maternal first heart sound, white Gaussian noise, environmental noises, and limited duration impulses considering as sensor noises. The epoch lengths were set equal to 430 and 400 ms for the analysis of the PCG drawn from the two

1 https://physionet.org/physiobank/database/simfpcgdb/. databases, respectively, using in both databases a window length of 50 samples.

#### **General Evaluation Indices**

The effectiveness of the WT–FD technique was tested *via* three general evaluation indices. The first *Q<sup>P</sup>* index calculates the efficiency of the algorithm in the correct detection of the S1 and S2 fetal heartbeat and its performance in the detection of locations that are not related with existing sounds. The *Q<sup>P</sup>* index is defined as follows:

$$Q\_P = 100\sqrt{\frac{\text{Sc}}{\text{S}\_O}} \frac{\text{S}\_C}{\text{S}\_P},\tag{19}$$

where *S<sup>O</sup>* is the number of sounds that every record contains, *S<sup>P</sup>* is the number of sounds that the proposed algorithm detects, and *S<sup>C</sup>* is the number of the *S<sup>P</sup>* sounds matching the *S<sup>O</sup>* sounds. Since the signals are simulated, the location of the existing sounds is specific, i.e., the *S<sup>O</sup>* number. For the *S<sup>C</sup>* number, the fHS was assumed to have been correctly detected when the estimated peaks lied in the intervals [*S*(2*i* + 1):*S*(2*i* + 2)], i.e., the start and end of each existing heart sound.

Furthermore, the second *D<sup>R</sup>* index indicates the percentage of the sounds that the WT–FD algorithm correctly detects out of the total number of sounds that it detects. The *D<sup>R</sup>* index is defined as follows:

$$D\_R = \frac{S\_C}{S\_P} 100.\tag{20}$$

Conclusively, the third *S<sup>F</sup>* index indicates the percentage of the sounds that the algorithm detects correctly out of the real fHS that every record contains. The *S<sup>F</sup>* index is defined as follows:

$$\text{S}\_F = \frac{\text{S}\_C}{\text{S}\_O} \text{100.} \tag{21}$$

The above three indices were calculated for the evaluation of the testing WT–FD method, and the results are presented in Section "Results".

#### **RESULTS**

As mentioned in Section "General Evaluation Indices," the WT– FD technique was tested on two simulated databases of fPCG signals. The results of this assessment are presented in **Tables 1** and **2**




**TABLE 2** | Performance of the wavelet transform–fractal dimension filter for signals of PhysioBank.

where the *QP*, *DR*, and *S<sup>F</sup>* indices are tabulated, providing a means for the evaluation of the performance of the WT–FD algorithm for the detection of the fetal S1 and S2 heartbeat and also each fHS separately.

In particular, **Table 1** presents the cases of 12 simulated fPCG signals created for the present study and consists 4 different fHRs and 3 different SNR values (white Gaussian noise). From **Table 1** it is clear that the WT–FD method is efficient for multiple conditions. The *Q<sup>P</sup>* index indicates that in all cases of fHR the WT–FD correctly predicts almost all the observed sounds in different SNR values up to 3 dB. Specifically, in cases of normal fHR (140 bpm), the algorithm has mean performance 100%. In cases of tachycardia (180 bpm), bradycardia (110 bpm), and arrhythmia, the efficiency of the method is slightly lower although it is sufficiently effective in the detection of the S1 beat locations.

Moreover, **Table 2** presents the cases of seven simulated fPCG signals from PhysioBank with different SNR values. Results for the cases of normal fHR with a range of SNR noise lying in [*−*26.3, *−*6.6 dB] demonstrate that the WT–FD algorithm segments and detects almost all the observed heart sounds and has a mean accuracy 81%. However, it is clear that the lower the SNR value, the harder it is for the WT–FD to segment and select the correct S2 fetal heartbeat. Very low SNR (less than *−*22.1 dB) makes the S2 sound difficult to distinguish from the noise. The *D<sup>R</sup>* index declares that, despite the fact that the algorithm misses a few heartbeats, it does not detect false locations. Most of the detected sounds are assigned to real fetal heartbeat locations. Furthermore, it is notable that all the detected S1 beat locations refer to real sounds.

**Figure 3** shows the efficiency of the WT–FD technique to recognize the fHS in signals with unexpected noise presence. **Figure 3A** corresponds to the *X*[*n*] unprocessed signal, and **Figure 3B** corresponds to the *X*REC[*n*] segmented reconstructed signal. In *X*[*n*] signal it is obvious that there is a noisy segment, which is marked with an arrow, that masked the S2 heart sound. In the *X*REC[*n*] signal it is clear that the WT–FD successfully extracts the sound.

The proposed WT–FD approach was also tested in real fPCG signals from a small pilot study, involving recordings from three pregnant women. The fPCG signals were recorded using vibration sensors (cost \$1 each) embedded in high definition 3Dprinted plastic harnesses. Each harness holds a ceramic piezo vibration sensor (35 mm diameter) on the maternal abdomen with rubber-made cushion to minimize the shear noise. The 3Dprinted harness is designed with precise parameters that rigidly

case with unexpected robust noise. **(B)** *X*REC[*n*] corresponds to the normalized treated signal without the overlap of noise. The arrows indicate the location of the S2 sound that the algorithm efficiently reveals.

mount the piezo sensor. Each sensor picks fPCG signals through a coaxial cable having very high insulating resistance. Power lab data acquisition system by AD instrument<sup>2</sup> was used to record the abdominal phonograms at a sampling frequency of *f<sup>s</sup>* = 1,000 Hz.

A characteristic example of one channel fPCG recording (time section of 3 s) with maternal heart rate of 96 bpm and fHR of 145 bpm is shown in **Figure 4A**. From the latter, it is clear that the fPCG signal is modulated by noise from various sources, and the most intense interferences are the mother's respiratory and heart sounds. **Figure 4B** shows the fourth level of the estimated WT coefficients from the eight level WT decomposition. The WT–FD method selects these WT coefficients that include information regarding the signal of interest, i.e., fHSs, based on the criterion (Eq. 9). For the real fPCG data processing, the constant of the

<sup>2</sup> http://www.adinstruments.com.

result of the de-noised fHS signal after the final WT–FD analysis with S1 and S2 denoting the first and second fHS, respectively.

criterion was set at 0.001, leaving out the first three decomposition levels and including only those with embedded fHSs. Finally, **Figure 4C** depicts the estimated fHS signal, i.e., the detected S1 and S2 fHSs, marked with (S1) and (S2), respectively, as the final output from the proposed WT–FD method. Note that, in some cases [**Figure 4C** around (0.5–1.5 s)], three S2 fHSs were missed by the WT–FD filter due to their lower intensity, compared to the neighboring S1 ones and the local background noise. Nevertheless, when comparing the original recording of **Figure 4A** with the outputted fHS signal from the WT–FD approach in **Figure 4C**, a clear contribution to the enhancement of the fHS signal from its original recording is evident.

It should be noted that the fHSs are not perfectly periodic due to the heart rate variability. It can be seen that, despite the noisy signal, the WT–FD method successfully identifies the fHSs and their time location and duration, giving the physicians the means

#### **REFERENCES**


to estimate the fHSs and the fHR. The wavelet morphology of the sounds could vary with different pathophysiological conditions. This is of great importance when the fHSs are continuously recorded for long-term analysis.

#### **CONCLUDING REMARKS**

The fPCG signals are of low amplitude and loaded with heavy noise. The sources of the noise, i.e., maternal sounds, fetal movement, sound produced by the transducer, and other, are overlapping the main fHS. The literature in the area of fetal auscultation is not strict about the intensity of background noise and the intensities of S1 and S2 heartbeat, because of the different auscultation devices but also due to the different gestation age. Nevertheless, it is possible to argue that the amplitude of the stationary background noise did not fully overlap the fHS, and that the SNR values that have been tested in the present study were sufficient samples of heavily loaded signals. However, as it was shown by the testing results, the WT–FD scheme is quite satisfactory in the analysis of the fHS. This first approach of the research in fPCG signals reveals sufficient information, which indicates that this technique can be a promising fHS segmentation tool. Furthermore, there are perspectives for low-cost and continuous recordings in homecare setups and diagnosis of conditions related, for example, to fetus maturation or specific abnormalities.

Future work will focus upon the extension of WT–FD to real recorded signals for a better review of fetal functionality and the fetal heart cycle. Moreover, multichannel recordings could be considered, taking into account the spatial orientation of the fetus and the proximity to the mother's heart sound noise. As phonocardiography has been an important field in the research area related to the fetus for some time, efficient characterization of fetal heartbeat could contribute to the automated determination of fetus parameters. In this vein, the determination of multiple fetus health data may reveal new aspects, which could improve the safety of pregnancies.

#### **AUTHOR CONTRIBUTIONS**

Producing the code used in the paper: EK and LH. Analyzing the signals and drawing conclusions as well as discussing the structure of the paper and writing up the final version: EK, LH, IC, and AK.

#### **ACKNOWLEDGMENT**

The real fPCG data were obtained from a research project funded by Abu Dhabi Education Council (ADEC).


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer, WT, and handling editor declared their shared affiliation, and the handling editor states that the process nevertheless met the standards of a fair and objective review.

*Copyright © 2017 Koutsiana, Hadjileontiadis, Chouvarda and Khandoker. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Efficient Fetal-Maternal ECG Signal Separation from Two Channel Maternal Abdominal ECG via Diffusion-Based Channel Selection

Ruilin Li <sup>1</sup> , Martin G. Frasch<sup>2</sup> \* and Hau-Tieng Wu1, 3 \*

<sup>1</sup> Department of Mathematics, University of Toronto, Toronto, ON, Canada, <sup>2</sup> Department of Obstetrics and Gynecology, University of Washington, Seattle, USA, <sup>3</sup> Mathematics Division, National Center for Theoretical Sciences, Taipei, Taiwan

There is a need for affordable, widely deployable maternal-fetal ECG monitors to improve maternal and fetal health during pregnancy and delivery. Based on the diffusion-based channel selection, here we present the mathematical formalism and clinical validation of an algorithm capable of accurate separation of maternal and fetal ECG from a two channel signal acquired over maternal abdomen. The proposed algorithm is the first algorithm, to the best of the authors' knowledge, focusing on the fetal ECG analysis based on two channel maternal abdominal ECG signal, and we apply it to two publicly available databases, the PhysioNet non-invasive fECG database (adfecgdb) and the 2013 PhysioNet/Computing in Cardiology Challenge (CinC2013), to validate the algorithm. The state-of-the-art results are achieved when compared with other available algorithms. Particularly, the F<sup>1</sup> score for the R peak detection achieves 99.3% for the adfecgdb and 87.93% for the CinC2013, and the mean absolute error for the estimated R peak locations is 4.53 ms for the adfecgdb and 6.21 ms for the CinC2013. The method has the potential to be applied to other fetal cardiogenic signals, including cardiac doppler signals.

Keywords: de-shape short time Fourier transform, fetal electrocardiogram, maternal abdominal electrocardiogram, nonlocal median, diffusion maps

# 1. INTRODUCTION

Fetal electrocardiogram (ECG) and the fetal heart rate (HR) provide enormous information about fetal health. For example, the fetal distress monitoring (Jenkins, 1989) or the potential risk for fetal hypoxia detection and alert by the ST analysis monitor (Belfort et al., 2015). Moreover, from clinical studies and animal models, evidence is accumulating that perinatal brain injury originates in utero, yet no means exist to detect its onset early, reliably and with simple, widely accessible means (Anblagan et al., 2016). A harbinger of brain injury is the fetal inflammatory response (Hagberg et al., 2015). There is an urgent need for early antenatal detection of fetal inflammatory response to prevent or at least mitigate the developing perinatal brain injury. In adults and neonates, complex mathematical features of heart rate fluctuations have proven promising as early diagnostic tools (Bravi et al., 2013; Fairchild et al., 2014). For the fetal monitoring, our team addressed the challenge by developing a series of biomarkers relying on non-invasively obtainable fetal HR. Our fetal inflammatory index tracks inflammation along with the fetal plasma IL-6 temporal profile in a fetal sheep model of subclinical chorioamnionitis (Durosier et al., 2015). We also derived a set of fetal HR features that is specific to brain or gut inflammation (Liu et al., 2016). Such systemic and organspecific tracking of inflammation via fetal HR is possible due to the brain-innate immune system

#### Edited by:

Ahsan H. Khandoker, Khalifa University, United Arab Emirates

#### Reviewed by:

Joachim Behar, Technion - Israel Institute of Technology, Israel Paolo Melillo, Università degli Studi della Campania "Luigi Vanvitelli" Caserta, Italy

#### \*Correspondence:

Martin G. Frasch mfrasch@uw.edu Hau-Tieng Wu hauwu@math.toronto.edu

#### Specialty section:

This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology

Received: 04 February 2017 Accepted: 18 April 2017 Published: 16 May 2017

#### Citation:

Li R, Frasch MG and Wu H-T (2017) Efficient Fetal-Maternal ECG Signal Separation from Two Channel Maternal Abdominal ECG via Diffusion-Based Channel Selection. Front. Physiol. 8:277. doi: 10.3389/fphys.2017.00277 communication reflected in the fetal HR fluctuations, commonly referred to as the cholinergic anti-inflammatory pathway (Fairchild et al., 2011; Olofsson et al., 2012; Garzoni et al., 2013).

In spite of its broad usefulness in the fetal health, it is fair to state that in the fetal HR monitoring realm, the technological progress has been coming more gradually. This has been not due to the plethora of studies attempting and testing various approaches, but, rather, due to the intrinsic limitations of the currently used fetal HR monitoring technology. This technology is outdated, as it deploys the traditionally set low sampling rate of heart rate or ECG signal. In animal model and human cohorts, we showed that such sampling rate is bound to miss the faster temporal fluctuations of vagal modulations of fetal HR variability and leads to inaccuracies in detection of early fetal acidemia (Durosier et al., 2014; Li et al., 2015). A sampling rate of the ECG signal around 1,000 Hz is required to capture these vagal influences and this is the commonly used sampling rate for the postnatal studies and our above-cited studies on the fetal inflammatory index.

Postnatal clinical studies are typically based on multi-lead ECG recordings which, even in newborns, and certainly in adults, poses no technical challenge to attach and record from. In fetuses, however, this is not the case. Since the fetal cardiac electric field strength is order of magnitude weaker than maternal ECG's, and the lack of clinical motivation in higher quality fetal HR data, little development had been done to focus on fetal ECG (fECG) signal in the clinical monitoring until today, except the Dopplerbased fetal HR extraction techniques that dominate the market. The Doppler-based fetal HR extraction techniques, however, suffer from low fetal HR sampling rates, largely due to the autocorrelation algorithms deployed in the devices (Durosier et al., 2014). Transabdominal ECG (aECG) machines overcome this limitation by capturing the actual cardiac electric field and have returned to the market during the last decade. However, their arrival has been slower than we would have hoped. Perhaps this is in part due to the general acceptance speed of new technology in medicine (related to regulatory and safety testing as well as the specific cultures), due to the high cost for each device to upgrade a hospital's delivery unit, or, more likely, the technical limitation of the fetal ECG extraction from the aECG signals. To make the technology of high quality and low-cost fetal ECG widely accessible, we need algorithms for fetal ECG extraction from easily deployable aECG devices. We refer the readers to the up-to-date guidance (Behar et al., 2016) for more information about the non-invasive fetal ECG extraction and analysis.

The current study addresses this challenge by proposing an algorithm capable of working with only two composite (maternal and fetal) aECG channels to derive the fetal signal from it. It is based on the currently developed single-lead fECG algorithm based on the modern time-frequency analysis and manifold learning technique (Su and Wu, 2017) and a novel proposed diffusion-based channel selection criteria. All the proposed methods have rigorous mathematical backups, and numerically they can be efficiently implemented to handle long signal. We call the proposed algorithm SAVER, which stands for **S**mart **A**dapti**V**e **E**cg **R**ecognition. To validate SAVER, we report the analysis results of two publicly available databases, and compare the algorithm with other available algorithms in the literature.

The paper is organized in the following way. In Section 2, we detail our proposed algorithm, describe the algorithms we will compare, and describe the databases we validate the algorithm. The results are shown in Section 3, and the discussions with the future works are provided in Section 4. The paper closes with the conclusion shown in Section 5. The necessary theoretical background is provided in SAVER Section SI.1 in appendix, particularly the diffusion-based channel selection criteria. We refer the readers to Su and Wu (2017) for the details of the deshape short time Fourier transform (dsSTFT), beat tracking and the nonlocal median.

# 2. METHODS

## 2.1. Two-Lead fECG Algorithm – **SAVER**

We now describe the proposed two-channel fECG algorithm, which the authors coined as SAVER. The overall algorithm is illustrated in **Figure 1**.

Denote two simultaneously recorded aECG signals as **x**0, **y**<sup>0</sup> ∈ R <sup>N</sup> with the sampling rate ξ0Hz over the interval from the 0-th second to the N/ξ0-th second. If the signal is sampled more slowly than 1,000 Hz, to enhance the R peak detection and the the nonlocal median (Su and Wu, 2017), the signal is upsampled to 1,000 Hz (Laguna and Sörnmo, 2000). We use the same notations to denote the upsampled signal.

#### **Step 0: pre-processing**

To suppress the noise, the signal is low-pass filtered below 100 Hz. Then, subtract the estimated trend from **x**0, **y**0, where the trends are estimated using median filter with window length LMF > 0 s. If needed, the power-line interference is suppressed by two notch-filters at 50 and 60 Hz, since the origin of the tested database in this paper is unknown (if the resource of the database is known, the notch-filter will be designed according to the power system of that region). Denote the pre-processed signal as **x** and **y**. Take a discrete finite subset I ⊂ (−1, 1]. Define **z**<sup>θ</sup> = θ**x** + √ 1 − θ <sup>2</sup>**y**, where θ ∈ I; that is, **z**<sup>θ</sup> is a linear combination of two aECG signals. This linear combination could be viewed as a generalization of the augmentation technique considered in Andreotti et al. (2014, Section 2.3.3).

#### **Step 1: maternal ECG estimation**

We iterate the dsSTFT and nonlocal median algorithms proposed in Su and Wu (2017) to decompose the maECG from each linear combination in {**z**<sup>θ</sup> }θ∈I. The algorithm is summarized below. For each θ we run the following three sub-steps.


signal, apply the beat tracking algorithm (Su and Wu, 2017, Section 3.1.3) to **z**θ to compute the locations of maternal Rpeaks. Denote the timestamps of estimated maternal R peaks as **r** m <sup>θ</sup> = (r m <sup>θ</sup>,1, . . . ,<sup>r</sup> m θ,kθ,<sup>m</sup> ), where kθ,<sup>m</sup> ∈ N is the number of estimated maternal R peaks.

3. (step 1-3) Adjust the estimated maternal R-peak locations by searching the maximum of **z**θ over a small window around **r** m θ . We use the same notation **r** m θ to denote the adjusted estimated maternal R-peak locations. Apply the nonlocal median (Su and Wu, 2017, Section 3.1.4) to estimate the maECG in **z**θ based on the estimated R-peak locations **r** m θ . Denote the estimated maECG as **z**˜θ,m.

#### **Step 2: channel selection**

For each linear combination in {**z**<sup>θ</sup> }θ∈I, with the estimated maECG, we obtain a rough fECG by a simple subtraction:

$$
\tilde{\mathbf{z}}\_{\theta,\underline{f}} \colon= \mathbf{z}\_{\theta} - \tilde{\mathbf{z}}\_{\theta,m}.\tag{2.1}
$$

Denote {˜**z**θ,<sup>f</sup> }θ∈<sup>I</sup> to be the set of rough fECG signals estimated from Step 1. We apply the lag map and the diffusion map (DM) to each rough fECG in {˜**z**θ,<sup>f</sup> }θ∈<sup>I</sup> and select the optimal linear combination by the following procedure. See Section SI.1 in the Appendix for the theoretical background of this approach.

For each rough fECG, say **z**˜θ,<sup>f</sup> , we evaluate the signal quality index (SQI) for the channel selection purpose in the following way. Apply the L-step lag map to embed the interval [2, TCS + 2] seconds of **z**˜θ,<sup>f</sup> into R L , where TCS > 0 is chosen by the user and 2 is chosen to avoid the boundary effect associated with the window in the dsSTFT approach. Here TCS is chosen to be short enough to guarantee the computational efficiency and to avoid the possibility non-stationarity inherited in the fECG signal, and long enough to capture the periodicity of the fECG. Denote the embedded point cloud as <sup>X</sup>θ,<sup>f</sup> ⊂ R L . Apply the 1-normalization DM to Xθ,<sup>f</sup> , where the bandwidth of the kernel is chosen in the following way suggested in Keller et al. (2010). We first set ǫ<sup>0</sup> to be the smallest value such that each data point has at least one neighbor within the distance ǫ0. Then we set the bandwidth to be 2ǫ0. Denote φθ,1 be the first nontrivial eigenvector of the corresponding graph Laplacian. Compute the power spectrum of φθ,1, denoted as |φˆ <sup>θ</sup>,1| 2 . Denote ξθ,1, ξθ,2, . . . , ξθ,nCS > 0, where nCS ∈ N is the number of peaks chosen by the user, to be the frequencies associated with the highest nCS peaks in |φˆ <sup>θ</sup>,1| 2 . Fix LCS > 0 and denote J<sup>θ</sup> : = ∪nCS j=1 [ξθ,<sup>i</sup> − LCS, ξθ,<sup>i</sup> + LCS]. The SQI for the channel selection purpose is thus defined as

$$\mathcal{S}\_{\theta} = \frac{\int\_{[0,\xi\_0/4)\cap\mathcal{J}\_{\theta}} |\hat{\phi}\_{\theta,1}(\xi)|^2 d\xi}{\int\_{[0,\xi\_0/2)\backslash\mathcal{J}\_{\theta}} |\hat{\phi}\_{\theta,1}(\xi)|^2 d\xi}. \tag{2.2}$$

Under the assumption that the better the quality of the rough fECG is, the closer the embedded point cloud is to the onedimensional circle, we know that the higher the SQI, the better the rough fECG is. More precisely, if the embedded point cloud is close to the one-dimensional circle, the first non-trivial eigenvector should behave like an oscillatory function. With the designed SQI, we could choose the optimal rough fECG as the one with the highest SQI. Denote **z**˜ ∗ f to be the optimal rough fECG with the highest signal quality index we can obtain from the given two channels.

#### **Step 3: fetal R peaks estimation**

With the rough fECG **z**˜ ∗ f obtained from the optimal linear combination, we finish the algorithm by estimating the fetal R peaks and fECG by again applying the dsSTFT and the nonlocal median algorithm. This part of the algorithm is essentially the same as that for the maternal ECG estimation, and we repeat the three sub-steps below for the sake of completeness.


Remark 2.1. We mention that by applying the nonlocal median again based on **r** f , we could denoise the optimal rough fECG waveform **z**˜ ∗ f and obtain a clean fetal waveform. However, since the result is similar to that shown in Su and Wu (2017), and the focus of this paper is the fetal R peak detection, we skip the details of the fECG reconstruction in this study, and leave the fetal waveform reconstruction in the future work.

#### 2.2. Comparison with Benchmark Algorithms

There have been several algorithms proposed in the field suitable for analyzing fECG from multiple channel aECG signals. Note that the two-channel aECG signals fall in the category of the blind source separation (BSS) (De Lathauwer et al., 2000; Akhbari et al., 2013; Di Maria et al., 2014; Varanini et al., 2014) and its variations (Sameni et al., 2008; Haghpanahi and Borkholder, 2013; Akbari et al., 2015). It is well known that usually we need more than 4 channels to have a reasonable result (Andreotti et al., 2016). Due to the stationarity assumption of the ICA, the input signal should be truncated to be short enough, like 30 s long. An important step in the BSS approach is channel selection, which is critical to identify the decomposed channel that contains the maternal or fetal ECG. Although we only have two channels, for the comparison purpose, we still show the results of the BSS approaches, including the joint approximation diagonalization of eigen-matrices (JADE) for the independent component analysis (ICA) and the principal component analysis (PCA). Since there are only two decomposed signals, we do not carry out the channel selection algorithms proposed in, for example, Andreotti et al. (2014); instead, we take the ground truth annotation to select the optimal channel that is more likely to be the fECG, and report the detected R peaks from this detected channel. Note that we do not take the ground truth annotation into account in any other algorithms considered in this paper except this BSS approach, due to the limited number of channels. We apply the publicly available codes provided in http://www.fecgsyn.com, and call the PCA method BSSPCA and the ICA method BSSICA following the terminologies suggested in Andreotti et al. (2016).

Another set of algorithms allow us to take only single mECG signal, but need to simultaneously acquire the maternal thoraciclead ECG signal (tECG). Examples include adaptive methods (AM) based on the least mean square (LMS) (Widrow et al., 1975) or the recursive least square (RLS) (Behar et al., 2014a) and its variations, like the echo state neural network (ESN) (Behar et al., 2014a), blind adaptive filtering (Graupe et al., 2008), extended Kalman filter (EKF) (Sameni, 2008; Niknazar et al., 2013; Andreotti et al., 2014), etc. In these algorithms, the maternal thoracic ECG signal (mtECG) is needed and is viewed as the reference channel. The mtECG contains the maternal cardiac activity information that we want to remove from the aECG. Based on the assumption that the mtECG and the maternal cardiac activity in the aECG are linearly related, the LMS or RLS helps to extract the fECG from the aECG by removing the maternal cardiac activity in the aECG. If the relationship between the tECG and the maternal cardiac activity in the aECG is nonlinear, then ESN could help. However, it is not always the case that we could get the mtECG, particularly in our setup, so these algorithms could not be directly applied for our purpose. Since it has been shown in Su and Wu (2017) that by combining the dsSTFT and nonlocal median, we are able to estimate the maECG signal accurately, we could thus view the estimated maECG signal as the reference. This consideration can also be found in, for example, Rodrigues (2014). For the LMS or ESN, we take the publicly available code from http://www.fecgsyn.com, and call the LMS method AMLMS and the ESN method AMESN following the terminologies suggested in Andreotti et al. (2016). We thus consider the following combinations of the proposed two channel fECG algorithm and the LMS or ESN. Precisely, in our proposed algorithm, we replace the direct subtraction (2.1) in Step 2 by the LMS or ESN, by taking the estimated maECG as the reference channel to get the rough fECG. We call the combined algorithm ds-AMLMS or ds-AMESN. Note that under the assumption that the nonlocal median does a good job to recover the maECG, the reference channel should be the same as, or linearly related to, the maternal cardiac activity in the aECG, so the LMS could be applied. The same idea could be applied to other algorithms, like RLS, but to keep the discussion simple, we focus on the above-mentioned two typical algorithms, LMS and ESN.

We could also consider the EKF algorithm. In the EKF algorithm, the information of the maternal R peak location is needed to cancel the maternal cardiac activity. Again, since it has been shown that by combining the dsSTFT and nonlocal median, we are able to estimate the maternal R peaks location accurately (Su and Wu, 2017), we could use the estimated maternal R peaks as the input to the EKF algorithm, and replace the direct subtraction (2.1) in Step 2 by the EKF. For the EKF, we take the publicly available code from http://www.fecgsyn.com, and call the EKF method TSEKF following the terminologies suggested in Andreotti et al. (2016). The combined algorithm is called the ds-TSEKF.

To have a complete comparison, we also consider the template subtraction (TS) algorithm, which is suitable for the single lead mECG signal, and replace the direct subtraction (2.1) in Step 2 by the TS algorithm. In this work, the TS method we apply is the singular value decomposition approach proposed in Kanjilal et al. (1997) and nominated in Behar et al. (2014b) and Andreotti et al. (2016) as TSPCA. For TSPCA, we take the publicly available code from http://www.fecgsyn.com. We call the combined algorithm ds-TSPCA. Other TS methods could be combined in the same way and we do not report the results to simplify the discussion.

We follow the suggested optimized parameters accompanying the code without any modification; for example, the input signal to the AM algorithms, like LMS or ESN, is resampled to 250 Hz<sup>1</sup> , the input signal to the TSEKF and TSPCA algorithms is resampled to 1000Hz, and we do not change the suggested initialization of the TSEKF code.

#### 2.3. Materials

We validate the proposed two-channel algorithm on two publicly available databases of aECG signals.

The first database is the PhysioNet non-invasive fECG database (adfecgdb), where the aECG signals with the annotation provided by experts are publicly available https://www.physionet. org/physiobank/database/adfecgdb/ (Goldberger et al., 2000; Kotas et al., 2010). There are five pregnant women between 38 and 40 weeks of pregnancy in this database. Each has 4 aECG channels and one direct fECG signal recorded from the Komporel system (ITAM Institute, Zabrze, Poland)<sup>2</sup> . The four abdominal leads are placed around the navel, a reference lead is placed above the pubic symphysis, and a common mode reference electrode with active-ground signal is placed on the left leg. See **Figure 2** for an illustration of the leads placement. The signal lasts for 5 min and is sampled at a fixed rate 1,000 Hz with the 16 bit resolution. The R peak annotation is determined from the direct fECG recorded from the fetal scalp lead.

The second database is the 2013 PhysioNet/Computing in Cardiology Challenge (https://physionet.org/challenge/2013/# data-sets), abbreviated as CinC2013. We focus on the set A composed of 75 recordings for an assessment of our proposed algorithm since it is the only one with the provided the R peak annotation with reference to a direct FECG signal, acquired from a fetal scalp electrode. Each recording includes four noninvasive mECG channels that were obtained from multiple sources using a variety of instrumentations with differing frequency response, resolution, and configurations. Although they are from different resources, all recordings are resampled at the sampling rate 1,000 Hz and last for 1 min. There is no publicly available information

<sup>1</sup>We mention that the publicly available code for AM algorithms is optimized for the 250 Hz signal, and it is likely that the results would be improved if the corresponding parameters had been optimized for 1,000 Hz.

<sup>2</sup>http://www.itam.zabrze.pl/developments-english-version-233/665-komporel.

about where the leads are placed on the maternal abdomen. Note that some recordings come from the adfecgdb database, but no detail is available publicly. More details about these two databases can be found on the website. We follow the suggestion in Andreotti et al. (2014) to disregard the recording a54 since it was discarded by the Challenge's organizers, and focus on the remaining 74 recordings.

#### 2.4. Evaluation Metrics

In the whole analysis, the R peak detection result is evaluated by beat-to-beat comparisons between the detected beats and the provided annotations. We follow the criterion in Guerrero-Martinez et al. (2006) and choose a matching window of 50 ms. Denote TP, FP, and FN to be true positive rate, false positive rate, and false negative rate, where TP means correctly detected peaks, FP means nonexistent peaks that were falsely detected, and FN means existing peaks that were not detected.

We report the sensitivity (SE) and the positive predictive value (PPV) defined as

$$SE \coloneqq \frac{TP}{TP + FN}, \quad PPV = \frac{TP}{TP + FP},\tag{2.3}$$

and the F<sup>1</sup> score, which is the harmonic mean of PPV and SE,

$$F\_1 \colon = \frac{2\,\text{TP}}{2\,\text{TP} + \text{FN} + \text{FP}}.\tag{2.4}$$

We also report the mean absolute error (MAE) of the estimated R peak locations. We follow the suggestion in Andreotti et al. (2016) to report the MAE only on true positive annotations to make the evaluation independent of the detection accuracy. Thus, the MAE is defined as

$$MAE \colon = \frac{1}{n\_{TP}} \sum\_{j=1}^{n\_{TP}} |r\_i^f - \tilde{r}\_i^f|, \tag{2.5}$$

where nTP is the number of true positive annotations, and r˜ f i and r f i are the temporal location of the i-th true positive reference Rpeak and temporal location of the i-th true positive detected R peak.

For each database, we will report two sets of statistics. First, for each subject, we record the best F<sup>1</sup> result among all pairs of available channels, denoted as F1(1) and report the mean and median of the F1(1) of all subjects, and the corresponding summary statistics of the MAE, denoted as MAE(1). To see how stable the algorithm is, we also record the median F<sup>1</sup> result among all pairs of available channels, called F1(0.5), and report the mean and median of the F1(0.5) of all subjects, as well as the corresponding summary statistics of the MAE, denoted as MAE(0.5). Second, to evaluate the lead placement issue, for each pair of available channels, we report the the mean and median of the F<sup>1</sup> of all subjects, and the corresponding summary statistics of the MAE. To avoid the boundary effect inevitable in the dsSTFT algorithm due to the window length, the first and last 2 s in every recording are not evaluated. The notation a ± b indicates the mean a with the standard deviation b.

#### 2.5. Parameters

For a fair comparison and the reproducibility purposes, here we summarize the parameters for SAVER. The parameters are fixed for all signals throughout the paper unless otherwise stated. For the linear combination of two channels, we fix I = {−1 + k/6} 12 k = 1 . The window length LMF of the median filter for the baseline wandering removal is chosen to be 0.1 s. For the dsSTFT, the beat tracking, and the nonlocal median, the parameters are set to be the same as those reported in Su and Wu (2017).

For the channel selection, we set the lag to L = 7 for the lag map; we choose the Gaussian kernel and α = 1 normalization for the DM; we choose TCS = 40, nCS = 6 and LCS = 0.1375 Hz for the adfecgdb database, and TCS = 10, nCS = 6 and LCS = 0.25 Hz for the CinC database. We mention that the above parameters are chosen in the ad-hoc fashion without any optimization pursue. Those parameters could be optimized based on the application field and the environment.

The algorithms are tested on MacBook Air (13-inch, Mid 2013) with Processor 1.3GHz Intel Core i5, Memory 4 GB1600MHz DDR3, Mac OS Sierra (Version 10.12.2), and Matlab R2015b without implementing the parallel computation.

#### 3. RESULTS

For the adfecgdb database, the direct fECG measurement was lost between 187 and 191 s and between 203 and 211 s in the r10 record, and these two segments were discarded in the evaluation. The evaluation results of our proposed algorithm for each combination of two channels out of four available channels of all subjects in the adfecgdb database are shown in **Table 1** for a clear comparison purpose. Except the combination of Channel 2 and Channel 3 in r01 and r08, all the other combinations have the F<sup>1</sup> consistently greater than 94%. For the MAE, the result is always smaller than 9ms except the combination of Channel 2 and Channel 3 in r08. **Table 2** shows the comparison



of the proposed method with other available algorithms. The F1(1) and F1(0.5) of all 6 pairs for each subject are recorded, and the summary statistics of all subjects are shown. It is clear that SAVER is consistently better than the other algorithms. The average running time is 141.55 s for SAVER, 194.44 s for the ds-AMLMS, 589.83 s for the ds-AMESN, 308.93 s for ds-TSEKF, 78.30 s for ds-TSPCA, 10.01 s for BSSICA, and 10.98 s for BSSPCA.

For the CinC2013 database, in **Table 3** we compare SAVER with the other available algorithms in the CinC2013 database. The F1(1) of all recordings of our method is 92.99 ± 16.0% and the corresponding MAE(1) is 5.38 ± 4.52 ms, which are both better than the other compared methods. The median F1(0.5) of all recordings of our method is 85.44±22.42% and the MAE(0.5) of our method is 6.54 ± 4.92 ms, which are both better than the best result determined by other methods. It should be noted that TABLE 2 | The summary statistics of different methods' performance, including F<sup>1</sup> and mean absolute error (MAE), evaluated in the adfecgdb database.


The F1(1) result from the six pairs of two channels is recorded for each subject, and the summary statistics of all subjects is reported in the top half rows; the F1(0.5) result from the six pairs of two channels is recorded for each subject, and the summary statistics of all subjects are reported from the bottom half rows. For the adaptive method (AM) part of ds-AMLMS or ds-AMESN, the TSEKF part of ds-TSEKF, and the TSPCA part of ds-TSPCA, we take the publicly available code from http://www.fecgsyn.com, and follow the suggested optimized parameters accompanying the code without any modification; for example, the input signal to the AM algorithms is resampled to 250 Hz, the input signal to the TSEKF and TSPCA algorithms is resampled to 1,000 Hz, and we do not change the suggested initialization of the TSEKF code. std, standard deviation; Q1, the first quartile; Q3, the third quartile.

the median of F1(0.5) over 6 pairs of our proposed algorithm is still as high as 96.32%, while other methods decline dramatically to less than 60%<sup>3</sup> . This result suggests the stability of the proposed

<sup>3</sup> It is suggested in Behar et al. (2014b) that the TSPCA algorithm performs better if we applied the 10 Hz high pass filter to remove the baseline wandering before applying the TSPCA algorithm. If we replace the median filter in Step 0 by the 10 Hz high pass filter to remove the baseline wandering, over all subjects except a54, we have F1(1) = 93.19 ± 14.19%, MAE(1)= 5.16 ± 3.85 ms, F1(0.5) = 74.79±27.74%, and MAE(0.5)= 9.26±6.34 ms. Note that while the results of F1(1) and MAE(1) are both better than all algorithms under comparison, the results of F1(0.5) and MAE(0.5) are worse. This discrepancy might be caused by the carried out optimization in Behar et al. (2014b). Since the baseline wandering algorithm is out of the scope of this paper, in all tables we only report the results with the median filter as mentioned in Step 0.


TABLE 3 | The summary statistics of different methods' performance, including F<sup>1</sup> and mean absolute error (MAE), evaluated in the CinC2013 database.

Method Mean Std Q<sup>1</sup> Median Q<sup>3</sup>

The subject a54 is removed from the datasets. The F1(1) result from the six pairs of two channels is recorded for each subject, and the summary statistics of all subjects is reported in the top half rows; the F1(0.5) result from the six pairs of two channels is recorded for each subject, and the summary statistics of all subjects are reported from the bottom half rows. For the adaptive method (AM) part of ds-AMLMS or ds-AMESN, the TSEKF part of ds-TSEKF, and the TSPCA part of ds-TSPCA, we take the publicly available code from http://www.fecgsyn.com, and follow the suggested optimized parameters accompanying the code without any modification; for example, the input signal to the AM algorithms is resampled to 250 Hz, the input signal to the TSEKF and TSPCA algorithms is resampled to 1,000 Hz, and we do not change the suggested initialization of the TSEKF code. std, standard deviation; Q1, the first quartile; Q3, the third quartile.

method<sup>4</sup> . The average running time is 20.29 s for SAVER, 27.26 s for the ds-AMLMS, 100.35 s for the ds-AMESN, 75.83 s for ds-TSEKF, 17.52 s for ds-TSPCA, 3.29 s for BSSICA, and 3.20 s for BSSPCA.

To further evaluate the influence of the lead placement, or to answer if we could design the best lead placement scheme for the proposed two-channel algorithm, we report the summary statistics of all pairs of two channels for the adfecgdb database in **Table 4** and the CinC2013 database in **Table 5**. It is interesting to see that for the adfecgdb database, except for the combination of channel 2 and channel 3, the mean F<sup>1</sup> accuracy is great than 97%. The outlier of the combination of channel 2 and channel 3 comes from the fact that the fECG is strong in case r08, which confuses the channel selection step. As a result, SAVER extracts the maternal ECG as the fECG, which leads to a wrong fECG estimation<sup>5</sup> . While determining the role of each component is a common issue for the fetal-maternal ECG separation algorithms and commonly we need more information to handle it, we leave this open problem for the future work.

Compared with the result of the adfecgdb database, the performance of SAVER in the CinC2013 database is not uniform cross different combinations of channels. Note that the lead placement scheme is unknown for the CinC2013 database, so it is not possible to conclude which pair of channels is the best. However, if we assume that the lead placement scheme for all recordings in the CinC2013 database is the same as the lead placement scheme shown in **Figure 2**, then the CinC2013 database results suggest that the best combination is channel 1 and channel 4; the F<sup>1</sup> has the mean of 87.93% with the standard deviation 22.64%, and the median 97.60% with the interquartile range 6.92%; the MAE has the mean of 6.21 ms with the standard deviation 6.03 ms, and the median 4.34 ms with the interquartile range 5.62 ms<sup>6</sup> . Another finding deserves a discussion is that unlike the adfecgdb database, we can see the discrepancy between the best F<sup>1</sup> out of the 6 pairs reported in **Table 3** and the average F<sup>1</sup> of each pair reported in **Table 5**. This might suggest that the lead system applied in the CinC2013 database is heterogenous across the recordings.

# 4. DISCUSSION

The encouraging results of SAVER indicate the possibility to design a "two-lead system" for the noninvasive, and long term fECG monitoring purpose. To the best of our knowledge, less is published about two aECG channels approach (for example, in Rodrigues, 2014, the considered algorithm can be applied to the two channel aECG), and our proposed method focuses on this direction. The main innovation of our approach, compared with other methods, is twofold. First, based on the geometry of the inherited oscillatory structure of the cardiac activity, the diffusion-based manifold learning technique is applied to do the

<sup>4</sup> It is suggested in Behar et al. (2014b, p. 1569) to remove six more recordings, a33, a38, a47, a52, a71, and a74, in addition to a54, because of some inaccurate reference annotations identified by the visual inspection of authors in Behar et al. (2014b). The F1(1) of all recordings of our method is 94.80 ± 13.17% and the MAE(1) is 5.04 ± 3.88 ms, and the F1(0.5) of all recordings of our method is 87.04 ± 21.27% and the MAE(0.5) of our method is 6.29 ± 4.56 ms.

<sup>5</sup> If we are allowed to use the physiological information that both the fetus and the mother are healthy so that the fetal IHR is on average higher than maternal IHR, then we could correct this confusion by swapping the fetal IHR and maternal IHR. This leads to the mean F<sup>1</sup> of the combination of channel 2 and channel 3 93.44% with the standard deviation 6.99% and the mean MAE 5.57 ms with the standard deviation 2.41 ms, and the results of other combinations unchanged.

<sup>6</sup> If we remove a33, a38, a47, a52, a54, a71, and a74 from the CinC2013 database (Behar et al., 2014b), for the combination of channel 1 and channel 4, the F<sup>1</sup> has the mean 89.81% with the standard deviation 20.84%, and the median becomes 98.41% with the interquartile range 5.10%; the MAE has the mean of 5.74 ms with the standard deviation 5.33 ms, and the median 4.20 ms with the interquartile range 5.48 ms.

TABLE 4 | The summary statistics of **SAVER**, including F<sup>1</sup> and mean absolute error (MAE), for six pairs of four available channels in the


std, standard deviation; Q1, the first quartile; Q3, the third quartile.

TABLE 5 | The summary statistics of **SAVER**, including F<sup>1</sup> and mean absolute error (MAE), for six pairs of four available channels in the CinC2013 database.


The subject a54 is removed from the datasets. std, standard deviation; Q1, the first quartile; Q3, the third quartile.

channel section. While other channel selection criteria mainly are based on the power spectral distribution, wave morphology entropy, root mean square error, etc, to find the clearest and most enhanced QRS complexes (Di Maria et al., 2014; Ghaffari et al., 2015), our approach is different since we carefully examine the nontrivial underlying geometric structure hosting the cardiac activity by the DM and look for the linear combination that is most like a simple closed curve. Second, we apply the modern time-frequency analysis technique, the dsSTFT, and the beat tracking algorithms detailed in Su and Wu (2017) to obtain an accurate R peak locations, and the nonlocal median, to better estimate the maternal ECG morphology and fetal ECG morphology. Compared with other available algorithms, we use more information hidden in the aECG, including decomposing the non-sinusoidal oscillatory pattern from the time-varying frequency, and the low dimensional parametrization of all possible cardiac oscillations. We mention that an important advantage of the approach in Su and Wu (2017) is the ability to separate mECG and fECG with temporal overlap by the nonlocal median. Furthermore, due to its nonlocal nature, it can directly handle a long signal without dividing it into small fragments. Notice that unlike the traditional AF-like methods, SAVER does not cancel the maternal ECG in one channel by designing a filter from another channel; instead, it directly cancels the maternal ECG in a single linear combination, as is mentioned in Step 1.

Our results deserve a discussion and comparison with the previous reported findings. For the adfecgdb database, our result is overall compatible with, or better than, the state-of-art result reported in the field. For example, if we choose the pair of channel 1 and channel 2, our result is better than the best channel result based on the continuous wavelet transform based single-channel algorithm (**Table 5**, Castillo et al., 2013). However, it is not a fair comparison since the algorithm used in Castillo et al. (2013) is a single-channel algorithm. On the other hand, if we compare with the methods based on ICA on four channels (**Table 1**, Poian et al., 2015), our result is compatible. The MAE, which is less reported in the literature, is as small as 10 ms, which indicates the potential of applying the SAVER to do the fetal heart rate variability (HRV) analysis.

For the CinC2013 database, our result is compatible, or better than, the reported results. At the first glance, it is not the case, since by the ICA-based algorithms (Andreotti et al., 2014; Behar et al., 2014b), the accuracy could be as high as have the mean F<sup>1</sup> = 96%, under the same setup that a detected R-peak was labeled as TP if within 50 ms of a reference R-peak. However, we mention that unlike SAVER, these algorithms are ICA-based and four channels are simultaneously used. Specifically, in Behar et al. (2014b, **Table 3**), among different combinations of different algorithms, the algorithm FUSE-SMOOTH achieved the best result – the mean F<sup>1</sup> over all recordings is 96%, after removing a33, a38, a47, a52, a54, a71, and a74; in Andreotti et al. (2014, **Table 1**), the augmentation, the ICA, the template adaptation or TSEKF, and other techniques are applied, and the result with the mean F<sup>1</sup> = 97.3% over all recordings with the standard deviation 0.108 is reported based on the template adaptation, after removing a54. Our proposed algorithm, on the other hand, outperforms the algorithm based on four channels and the BSSPCA, for example, Di Maria et al. (2014). In Di Maria et al. (2014, Section 3.2), the accuracy of the proposed algorithm in detecting the fetal heart beats gives the mean F<sup>1</sup> = 89.8% over all recordings, under the setup that a detected R-peak was labeled as TP if within 100 ms of a reference R-peak and removing 9 recordings, including a29, a38, a54, a56, a33, a47, a52, a71, and a74. Another novel method based on the channel selection over 4 channels followed by the sequential total variation denoising (Lee and Lee, 2016, **Table 5**) leads to the accuracy with F<sup>1</sup> = 89.9% and the MAE = 9.3 ms<sup>7</sup> under the setup that a detected R-peak was

<sup>7</sup> In a private communication, the authors confirmed that this F<sup>1</sup> is the "overall F1," which is evaluated by collecting all beats from all recordings, and evaluate the F<sup>1</sup> on all collected beats. If we follow the same procedure and remove a33, a38, a47, a52, a54, a71, and a74, the overall F<sup>1</sup> of SAVER for the combination of channel 1 and channel 4 is 89.77% and the MAE is 4.91 ms.

labeled as TP if within 50 ms of a reference R-peak and removing a33, a38, a47, a52, a54, a71, and a74. We emphasize that while our algorithm does not outperform some of the above-mentioned algorithms, based on two channels, SAVER leads to the MAE as small as 6.21 ms in channel 1 and channel 4 combination in the CinC2013 database, which again indicates the potential of applying the SAVER to do the fetal HRV analysis.

As discussed above, theoretically, the chance is low that the fetal cardiac axis orientation would be so much orthogonal to the 2-dim affine subspace spanned by the two leads that no fECG shape can be reconstructed. This is a big advantage compared with the single-lead system, as the chance that the fetal cardiac axis orientation is orthogonal to the 1-dim affine subspace spanned by the single lead is much higher. Thus, while there have been several successful algorithms for the one aECG channel, like (Castillo et al., 2013; Behar et al., 2014a; Su and Wu, 2017) and the citations inside, if the recorded one channel signal does not have fECG information, there is nothing the algorithm can do. From the practical viewpoint, since only two leads are needed, the corresponding hardware could be lighter and more deployable than the currently available four-lead or multiplelead systems. While it is certainly possible to generalize our algorithm to a three-lead or four-lead system (and the algorithm can be changed directly according to the setup), to have a better balance between the prediction accuracy, the hardware design, and practical purposes, we focus on the two-lead system in our research.

Despite of the above-mentioned benefits, there are several challenges we need to solve until this possible system is clinically usable. As is shown above, the performance of SAVER depends on how the two leads are put on the abdomen. The fECG situation is clearly different from the adult ECG system, like the widely applied 12 lead ECG system. Since fetus does move and rotate inside the uterus, the uterus differs from female to female, and the maternal body profile varies, we may not expect to have a two-lead system universal for all women. Therefore, for the practical purpose, particularly for the long term monitoring purpose and the future digital health, like the wearable biosensors (Li et al., 2017), it is important to ask if we could adaptively find the best lead placement scheme for different females. For the practical purpose, due to the inevitable non-stationary noise of different types, like the motion artifact and uterine contraction, an automatic system providing a SQI to alarm/warn the low quality of the lead system, and hence improve the overall fECG extraction quality, is urgently needed. We leave this important engineering problem to the future work. Another interesting question naturally raises from the current work is if we could generalize the current algorithm to study the twin dataset. Theoretically it is possible, if we take the fact that geometrically the twin will locate in different positions. We would expect to study this problem when the dataset is available.

From the algorithmic viewpoint, there are several directions we could improve the proposed two-channel fECG algorithm. The main ingredient in SAVER is the diffusion geometry. Since we have more than one aECG channel, we could consider modern diffusion-based manifold learning techniques to extract information common in two channels, like the alternating diffusion (Papyan and Talmon, 2016; Talmon and Wu, 2016; Lederman and Talmon, in press). The non-stationary nature of the fECG signal, which often presents itself as a time-varying frequency, might jeopardize the diffusion-based approach. We could consider to entangle the nontrivial timevarying frequency nature of the signal by further applying the modern nonlinear-type time-frequency analysis technique, like the synchrosqueezing transform or concentration of frequency and time (see Daubechies et al., 2016 and the citations inside). In this work, the parameters for the channel selection are chosen in the ad hoc fashion and are fixed across different algorithms for a fair comparison. For the practical purpose, we may optimize these parameters to improve the results. A systematic survey of this issue will be reported in the future work.

Another important algorithmic question left unanswered in this paper is how to improve the nonlocal median algorithm so that the reconstructed fECG could provide more accurate electrophysiological information about the heart, for example, the ECG morphology like the Q wave and ST-segment section information (Amer-Wåhlin et al., 2001). The main difficulty encountered in this problem is the lack of the "ground truth," and a careful design of the clinical trial to acquire a reliable ground truth for the morphological study of the fetal cardiac activity is needed. As important as this clinical information could be, we will focus on it as an independent research and report the result in the future work.

Last but not the least, the databases we tested are small and not specifically designed for our purpose. We thus need a large scale and well designed prospective study to confirm the result.

# 5. CONCLUSION

A novel two-channel fetal-maternal ECG signal separation algorithm, SAVER, is proposed. The potential of the proposed algorithm is supported by the positive validation results on two publicly available databases. The algorithm is both computationally efficient and is supported by the underlying rigorous mathematical model and theory. Its clinical applicability will be evaluated in the future work.

# AUTHOR CONTRIBUTIONS

RL and HW carried out the data analysis. MF and HW carried out the paper writeup.

# ACKNOWLEDGMENTS

HW research is partially supported by Sloan Research Fellow FR-2015-65363. The authors acknowledge the help of Jan Hamanishi to prepare the illustration in **Figure 2**. The authors thank the reviewers for their constructive comments to improve the manuscript.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys. 2017.00277/full#supplementary-material

### REFERENCES


associated with abnormal neuroimaging and outcomes in extremely low birth weight infants. J. Perinatol. 34, 375–379. doi: 10.1038/jp.2014.18


neonate? Pediatr. Crit. Care Med. 17, 165–176. doi: 10.1097/PCC.0000000000 000643


**Conflict of Interest Statement:** The authors filed a patent for this method.

Copyright © 2017 Li, Frasch and Wu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Evaluation of Abdominal Fetal Electrocardiography in Early Intrauterine Growth Restriction

Clarissa L. Velayo<sup>1</sup> \*, Kiyoe Funamoto<sup>2</sup> , Joyceline Noemi I. Silao<sup>3</sup> , Yoshitaka Kimura<sup>2</sup> and Kypros Nicolaides <sup>4</sup>

*<sup>1</sup> Department of Physiology, College of Medicine, University of the Philippines, Manila, Philippines, <sup>2</sup> Department of Obstetrics and Gynecology, Graduate School of Medicine, Tohoku University, Sendai, Japan, <sup>3</sup> Department of Obstetrics and Gynecology, Philippine General Hospital, University of the Philippines, Manila, Philippines, <sup>4</sup> Harris Birthright Research Centre, Kings College Hospital, London, United Kingdom*

#### Edited by:

*Joseph L. Greenstein, Johns Hopkins University, United States*

#### Reviewed by:

*Rathinaswamy Bhavanandhan Govindan, Children's National Medical Center, United States Miles J. De Blasio, Baker Heart and Diabetes Institute, Australia*

> \*Correspondence: *Clarissa L. Velayo clvelayo@up.edu.ph*

#### Specialty section:

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology*

> Received: *18 March 2017* Accepted: *09 June 2017* Published: *26 June 2017*

#### Citation:

*Velayo CL, Funamoto K, Silao JNI, Kimura Y and Nicolaides K (2017) Evaluation of Abdominal Fetal Electrocardiography in Early Intrauterine Growth Restriction. Front. Physiol. 8:437. doi: 10.3389/fphys.2017.00437* Objectives: This descriptive study was performed to evaluate the capability of a non-invasive transabdominal electrocardiographic system to extract clear fetal electrocardiographic (FECG) measurements from intrauterine growth restricted (IUGR) fetuses and to assess whether abdominal FECG parameters can be developed as markers for evaluating the fetal cardiac status in IUGR.

Methods: Transabdominal FECG was attempted in 20 controls and 15 IUGR singleton pregnancies at 20+0−33+6 weeks gestation. Standard ECG parameters were compared between the study groups and evaluated for their correlation. Accuracy for the prediction of IUGR by cut off values of the different FECG parameters was also determined.

Results: Clear P-QRST complexes were recognized in all cases. In the IUGR fetuses, the QT and QTc intervals were significantly prolonged (*p* = 0.017 and *p* = 0.002, respectively). There was no correlation between ECG parameters and Doppler or other indices to predict IUGR. The generation of cut off values for detecting IUGR showed increasing sensitivities but decreasing specificities with the prolongation of ECG parameters.

Conclusion: The study of fetal electrocardiophysiology is now feasible through a non-invasive transabdominal route. This study confirms the potential of FECG as a clinical screening tool to aid diagnosis and management of fetuses after key limitations are addressed. In the case of IUGR, both QT and QTc intervals were significantly prolonged and thus validate earlier study findings where both these parameters were found to be markers of diastolic dysfunction. This research is a useful prelude to a test of accuracy and Receiver Operating Characteristics (ROC) study.

Keywords: fetal electrocardiography, intrauterine growth restriction, prenatal screening, fetal cardiophysiology, fetal monitoring

**48**

# INTRODUCTION

Intrauterine growth restriction (IUGR) affects 3–10% of all pregnancies and has been associated with the development of adult cardiovascular diseases.(Barker et al., 1989; Bahtiyar and Copel, 2008) Recent studies suggest that cardiovascular changes may be present during the prenatal period with evidence of cardiac remodeling in the form of increased transverse diameters and more globular cardiac ventricles (Crispi et al., 2010). These corroborate with earlier descriptions of fetal hemodynamic adaptive changes in IUGR (Verburg et al., 2008) wherein the standard assessment of cardiac function utilizes various biochemical and ultrasound indices in relation to disease severity (Crispi et al., 2008). Until now, very few studies have successfully utilized electrocardiographic parameters to elucidate ongoing fetal cardiac dysfunction.

Our study is based on the intimate relationship between the electrical conduction system and the overall structure of the heart. Since IUGR induced cardiovascular remodeling affects both cardiac structure and function, and electrocardiographic parameters are established indicators of cardiac function as well, a thorough investigation of fetal electrocardiography may help us understand the pathophysiology behind fetal growth restriction.

In 1906, Cremer inadvertently discovered transabdominal fetal electrocardiographic (FECG) extraction (Skillern et al., 2005). Many physicians since then have affirmed its potential clinical value. However, non-homogenous tissue conduction, insulating effects of the vernix caseosa, poor signal noise reduction and its inability to pinpoint the exact cardiac structural changes prevented its acceptance for clinical use. The progress in the development of abdominal FECG systems then lagged most especially after the development of transcervical extraction techniques that addressed the issues raised with the transabdominal route. The transcervical route through fetal ECG ST segment analysis (STAN) is currently being used in intrapartal monitoring. The S-T interval on ECG represents myocardial repolarization. Its use in conjunction with cardiotocography was found to significantly reduce the number of fetal scalp blood sampling, operative vaginal deliveries and Cesarean section for the indication of fetal distress. European and North American trials (Vayssiere et al., 2007) showed good correlation between STAN and fetal cord blood pH. One study reported that a significant ST event for screening pH < or = 7.15 (21.9%) was with sensitivity 38% (41/108), specificity 83% (252/303), PPV 45% (41/92) and NPV 79% (252/319), and for pH < or = 7.05, it was (3.4%), 62.5% (10/16), 79% (313/395), 11% (10/92), and 98% (313/319), respectively. In another study, a negative predictive value of 95.2% was achieved (Devoe et al., 2006). Other studies have found FECG useful in: (a) evaluating atrioventricular block (AVB) in Anti-RO positive pregnancies, (b) monitoring drug therapy, or (c) monitoring fetal anemia (Gardiner et al., 2007).

Recent advances in the technology of abdominal FECG addressed previous limitations thus prompting the reintroduction of transabdominal FECG (Sato et al., 2007). A novel system developed in Sendai, Japan was used in this study. The basic advantages of this system over other techniques were: (1) it was non-invasive, lead electrodes were applied on the maternal abdomen using an adhesive pad; (2) it eliminated interference from maternal cardiac signals and fetal movement thus producing clear ECG recordings; and (3) it compared well with the gold standard, magnetocardiography, in terms of gestational age-dependent ECG parameters, such as the QRS, PQ and QT intervals (Horigome et al., 2000). A minimum recording time of 20 min was usually needed from which the fetal heart rate could be detected in real time and online in more than 70% of patients. In about 30%, 5–30–min off line calculations were required. For this study, all results were calculated off line to ensure good resolution of signals. Moreover, the inclusion of simultaneous Doppler studies increased the resolution of the FECG extraction process.

Fetal ECG (FECG) has a voltage of 5–20µV and signal amplitude of 1/50 of the maternal ECG (MECG) complex which averages 1,000µV. For this study, the method used was blind source separation with reference signals (BSSR) wherein the FECG was isolated from the surrounding noise by bouncing the signal through various transfer points. From a mathematical point of view, FECG is derived from the following linear equations:

ChA = a1 • FECG + a2 • MECG + noise1 ChB = b1 • FECG + b2 • MECG + noise2 ChC = c1 • FECG + c2 • MECG + noise3 where Ch = channel FECG = fetal electrocardiography MECG = maternal electrocardiography

Measured signals like ChA, ChB, and ChC are mixed signals that include FECG, MECG and other sources of noise (Kimura et al., 2012). In the presence of large noise, reference signals, such as continuous Doppler recordings or adaptation of a lead ECG waveform allows for shape mimicry of FECG and for detection by imitating signal timing. Moreover, our method differs from other systems because it is non-invasive and mobile. It is non-invasive since it merely necessitates placement of adhesive electrodes on the maternal abdomen. In contrast ST analysis (STAN) (Ross et al., 2004), requires rupturing the amniotic membranes and transcervically attaching needles into the fetal scalp. Mobility is achieved using a simple ensemble of portable equipment, laptop and signal extractor, unlike other immobile complex systems, such as in fetal magnetocardiography (FMCG) (Leuthold et al., 1999; Kähler et al., 2002; Grimm et al., 2003).

The cardiac conduction system and the structural morphology of the heart are intimately related thus making fetal electrocardiography a possible means to assess alterations in cardiac structure and performance. Cardiac overload or increases in peripheral resistance (McDade et al., 2001; Mari et al., 2007; Romo et al., 2009; Bolin et al., 2016) are known to occur in cases of fetal growth restriction. The objective of this study was to assess whether abdominal FECG parameters can be developed as markers for evaluating the fetal cardiac status in IUGR.

#### TABLE 1 | Study criteria.

#### INCLUSION CRITERIA EXCLUSION CRITERIA

	- A. Stage 0 includes fetuses with normal umbilical artery and middle cerebral artery Doppler indices.
	- B. Stage I includes fetuses with abnormal Doppler indices of the umbilical artery (increased resistance indices) or middle cerebral artery (cerebroplacental ratio of < 1).
	- C. Stage II includes fetuses with absent or reversed Doppler flow of the umbilical artery.
	- D. Stage III includes fetuses with absent or reversed Doppler flow of the a-wave of the ductus venosus.

## MATERIALS AND METHODS

In the Philippines and United Kingdom, all pregnant women who consulted or were admitted in the participating institutions, and were referred for fetal wellbeing studies, had met the inclusion criteria, and who did not have any of the exclusion criteria were included in the study (**Table 1**). The principal investigator obtained informed consent only after a thorough explanation and discussion of the study with the patient. Each patient then underwent a physical examination and interview prior to FECG recording and standard ultrasonographic and Doppler evaluation.

FECG was recorded using a novel extraction method from composite abdominal signals involving cancellation of the maternal ECG and blind source separation with a reference signal (BSSR). The method of abdominal FECG was as follows (Kimura et al., 2006): A standardized placement of 14 adhesive electrodes was utilized:12 placed on the maternal abdomen, inclusive of one reference electrode; one placed at the right maternal thoracic position; and one placed on the maternal back. Bipolarly recorded data of channels were sampled every 1 ms at 1 kHz with a 16-bit resolution which was subjected to 1–100 Hz band-pass filtering. Simultaneous fetal cardiotocography was performed to identify the opening and closing of fetal cardiac valves.

Each FECG recording lasted 20 min. All recordings were anonymized and sent to Sendai, Japan for fetal ECG extraction. For each recording, an averaged ECG waveform was used as a reference signal for the BSSR algorithm to seek fetal QRST components in the orthogonal signals from the first and second channel after the BSSR. This is followed by fast non-linear state space projection (FNSSP) within the time domain to precisely separate noise from the fetal ECG. Also, the simultaneous Doppler signals measured represented the opening and closing of aortic valves and were used to distinguish the end of the T wave. The fetal ECG standard parameters evaluated included the following: QT, RR, QRS, ST, PR, and PQ intervals, height of P wave, QTc, PR/RR and heart rate. Raw data was returned to the principal investigator for statistical analysis. The correlation of these standard fetal ECG parameters with other ultrasound [amniotic fluid deepest pocket (AF DP) or index (AFI), head circumference to abdominal circumference ratio (HC/AC), femoral length to abdominal circumference ratio TABLE 2 | Comparison of patients' demographics between CONTROL and IUGR groups.

changes 6. Severe maternal illness

The presence of any of the following: 1. Chromosomal anomalies 2. Congenital anomalies 3. Premature labor

4. Maternal medical contraindications to the use of electronic devices, such as pacemakers 5. Maternal mental disability, coma or sensorial


*p* < *0.05 alpha* \**(significant). SD, standard deviation.*

(FL/AC)] and Doppler markers [middle cerebral artery (MCA), umbilical artery (Umb A), uterine arteries (UA), and ductus venosus (DV)] were analyzed.

Data were collated and encoded in MSEXCEL 2013. Categorical data in the demographics, the maternal characteristics, and the diagnostic tests were described using frequency and percentages. Moreover, the continuous types of data were expressed in mean and standard deviation.

In comparing the mean and variances between control and FGR groups, the unpaired t-test was utilized to check the equality of variances as well as differences of means. Moreover, in testing associations among categorical types of data and groupings, the Chi Square test of independence and the 2 × 2 Fisher exact test were performed. Any associated p < 0.05 alpha were considered significant.

Furthermore, all the abdominal fetal electrocardiography parameters were tested for accuracy, such as sensitivity, specificity, negative, and positive predictive values. Any cut off values found to produce accuracy indices were presented.

The accuracy of computations was ensured through careful data processing and aided by IBMSPSS version 21 and NCSSPASS 2000 statistical software.

#### RESULTS

A total of 46 fetal ECG recordings were collected throughout the study period, 16 from the Philippines (RP) and 30 from the United Kingdom (UK). Out of the 46, only 20 Control and 15

June 2017 | Volume 8 | Article 437

IUGR fetuses were eligible after the final exclusion of patients with incomplete clinical data or poor ECG recordings.

Patients' demographics, maternal characteristics and clinical profiles were summarized in **Tables 2**–**4**. The mean age of patients was 30.86 ± 5.75 years (range: 17–40 years) where a majority (54%) of the study subjects were Asian (p = 0.018). All pregnancies (100%) were achieved through spontaneous conception in women with a mean gravidity score of 2.26 ± 1.52 and only 10 (28.6%) of these women reported a previous fetal loss. This evaluation revealed that there was a significant difference between the maternal weights of Control (C) and

TABLE 3 | Comparison of maternal characteristics between CONTROL and IUGR groups.


*p* < *0.05 alpha* \**(significant).*

*HT, height; WT, weight; BMI, body mass index, GA, gestational age.*


TABLE 4 | Comparison of clinical profiles between CONTROL and IUGR groups.

*p* < *0.05 alpha* \**(significant).*

*EFW, estimated fetal weight; HC, head circumference; AC, abdominal circumference; FL, femoral length.*

IUGR fetuses wherein mothers of the latter, on the average, weighed less (p = 0.034). Of all the indices of IUGR, HC/AC ratio alone was significant in predicting IUGR (p < 0.01). Moreover, for the IUGR sample population studied, there was an even distribution among the various stages of IUGR (11–17%) except in Stage III (3%).

Diagnostic Test Results (**Table 5**) showed significant differences between the C and IUGR groups. The most relevant changes were seen in the umbilical artery pulsatility index (Umb A PI), right and left uterine artery PI (UA PI) measurements and amniotic fluid deepest pool (AF DP) whose p-values were equal to 0.001, 0.001, < 0.01, and 0.040, respectively.

TABLE 5 | Comparison of diagnostic test results between CONTROL & IUGR groups.


*p* < *0.05 alpha* \**(significant).*

*PI, pulsatility index; EDF, end diastolic flow; MCA, middle cerebral artery; Vmax, maximum velocity; CP, cerebroplacental ratio; DV, ductus venosus; RT UT, right uterine artery; LT UT, left uterine artery.*

Clear P-QRST complexes were recognized in all cases (**Table 6**, **Figure 1**). In the IUGR fetuses, both the QT and QTc parameters were significantly prolonged (p = 0.017 and p = 0.002, respectively). There was no correlation between significant ECG parameters and Doppler or other indices to predict IUGR (**Tables 7, 8**). A detailed table matrix of cut off values in detecting IUGR per fetal ECG parameter was generated and showed an inverse correlation, increasing sensitivities but decreasing specificities, with the prolongation of ECG parameters (**Table 9**).

#### DISCUSSION

Intrauterine growth restriction (IUGR) is defined as the failure of a fetus to reach its growth potential (Mari and Hanif, 2008). It is

TABLE 6 | Comparison of fetal electrocardiographic parameters between CONTROL and IUGR groups.


*p* < *0.05 alpha* \**(significant).*

*QTc-corrected QT interval; HR, heart rate.*

distinguished from fetuses that are merely small for gestational age (SGA), weighing below the 10th percentile, by signs of chronic hypoxia. It is a condition caused by several factors, alone or in combination, such as maternal disease, malnutrition, drugs and toxins, uteroplacental insufficiency, fetal aneuploidy, fetal abnormality, genomic imprinting, and congenital infections. This can be further classified as LATE (TERM or ASYMMETRIC) IUGR, which includes fetuses with onset of growth restriction at a gestational age ≥34 weeks or EARLY (PRETERM or SYMMETRIC) IUGR, wherein growth restriction occurs at <34 weeks gestation. This study focused on early IUGR fetuses for two reasons: first, the non-invasive fetal ECG system was best calibrated for fetuses between 16 and 32 weeks age of gestation and second, it was during this exact period, <34 weeks gestation, that the detection and evaluation of IUGR was most critical to reduce perinatal morbidity and mortality secondary to iatrogenic or indicated premature delivery.

IUGR is brought about by a persistent state of low oxygen or chronic hypoxia from increased placental impedance. This leads to the shunting of blood in the peripheral circulation toward more vital organs or the "Brain Sparing Effect." Therefore, the fetal heart will undergo both physiologic and anatomic adaptation to meet the demands of increasing cardiac output. This constant cardiac remodeling eventually leads to both systolic and diastolic dysfunction where significant compensatory changes may be reflected in electrocardiographic as well as Doppler variations.

Overall, an increase in heart size due to free wall hypertrophy without ventricular dilatation causes progressive hemodynamic deterioration in IUGR fetuses (Bahtiyar and Copel, 2008). This is exemplified by a predictable pattern of worsening indices which progressively unfolds: umbilical artery pulsatility index increase, middle cerebral artery pulsatility index decrease, changes in right diastolic indices (right E/A, ductus venosus), changes in right systolic indices (right ventricular ejection force), then finally, both left diastolic and systolic cardiac dysfunction.

FIGURE 1 | Sample fetal ECG recording in an IUGR fetus. (A) Shows the actual fetal ECG of an IUGR fetus at 34 weeks gestation. (B) Shows the averaged waveforms of the same interval. (C) Shows the Doppler wave signal recorded simultaneously by an attached continues Doppler transducer. (D) Shows how to measure the T wave, wherein Ao, opening of aortic valve; Ac, closing of Aortic valve.

TABLE 7 | Correlation of fetal electrocardiographic parameters and IUGR indices in the CONTROL group.


*p* < *0.05 alpha* \**(significant).*

*AFI–amniotic fluid index; AF DP – amniotic fluid deepest pocket; Umb A–umbilical artery.*

TABLE 8 | Correlation of fetal electrocardiographic parameters and IUGR indices in the IUGR group.


*p* < *0.05 alpha* \**(significant).*

In clinical practice, Doppler indices in IUGR denote fetal cardiovascular status as well as the placental vascular conditions. Other standard measurements (AFI and the biometric parameters) on the other hand, are important signs of ongoing disease. The discovery of any significant correlation between ECG measurements and other IUGR markers, although not evident in this study, would allow a better understanding of the initiation and progression of cardiac structural change and its deterioration.

Our study findings show that fetal electrocardiography can distinguish between normal and IUGR fetuses as exhibited by the significant prolongation of QT and QTc intervals in the latter study group. Both these electrocardiographic parameters represent depolarization and repolarization of the ventricles and have been used as an indicator for deteriorating ventricular performance (Martin et al., 1994). These findings were consistent with earlier published data by Velayo et al evaluating abdominal fetal electrocardiography in congenital heart defects (CHD)


TABLE 9 | Test of accuracy for the prediction of IUGR in cut off values of the different fetal electrocardiographic parameters.

TABLE 9 | Continued


in singleton pregnancies wherein certain cardiac pathologies exhibited particular ECG changes (Velayo et al., 2011). For example, PR and QTc intervals were prolonged in fetuses with dilated cardiomyopathy (DCM). They concluded that FECG had a role in the analyses of fetal cardiac pathophysiology not only in congenital heart disease but also in acquired cardiac maladaptive geometry which is the case in IUGR (Saba et al., 2005). The latter also supports the results of studies on twin to twin transfusion syndrome where progression of disease and simultaneous fetal cardiac changes can be evaluated. The same research group in collaboration with a center in the United Kingdom published an evaluation of fetal cardiac performance in twin-to-twin transfusion syndrome wherein the cardiovascular status of the donor and recipient twins differed based on fetal Doppler studies (Velayo et al., 2012). Again, the QT interval and QTc were significantly prolonged this time in the donor fetus. This fact was interesting since compared to its enlarging recipient counterpart, the evolution of a donor twin has been likened to a fetus with IUGR.

Studies using FMCG (Leuthold et al., 1999; Kähler et al., 2002; Grimm et al., 2003; Bolin et al., 2016) have shown that in the normal fetus ECG parameters lengthen as gestational age increases. The prevailing theory for this is that prolongation is caused by the increase in size of the fetal heart. Moreover, in FMCG studies involving IUGR fetuses, P wave, PR and QRS were significantly shorter, making the abbreviation of parameters signs of cardiac abnormality (Grimm et al., 2003; Bolin et al., 2016). Our findings using an abdominal FECG system agree with the first conclusion, where the lengthening of intervals with age is dependent on the cardiac mass. But our results of prolonged ECG parameters in this study and in earlier reports by our group, are inconsistent with FMCG findings (Velayo et al., 2011, 2012). A possible explanation for the prolongation in parameters is that in IUGR fetuses, the cardiac size increases in relation to fetal size, part of the brain sparing effect.

Interestingly, recent reports involving other non-invasive abdominal FECG systems have shown promising results for the use of phase-rectified signal averaging measurements (PRSA) and fetal heart rate variability (FHRV) in IUGR fetuses (Stampalija et al., 2015a,b). PRSA signals are representative of the current maturation status of the fetal autonomic nervous system (ANS) where lower acceleration capacity (AC) and deceleration capacity (DC) signals have been observed in IUGR fetuses with or without brain sparing. One study found a significant association between

*(Continued)*

either AC and DC and the middle cerebral artery pulsatility index (PI; p = 0.01; p = 0.005) (Stampalija et al., 2015a). With FHRV, short term variability (STV) evaluation is known to reflect the status of the fetal oxygen supply, therefore indicating the fetal acid-base status. STV was seen to shorten in cases of IUGR (Stampalija et al., 2015b). All these parameters can identify the different behaviors of the heart rate and will give investigators a panoramic view of electrocardiopathology.

In future investigations, we propose the inclusion of parameters, such as the cardiac circumference (CC) to thoracic circumference (TC) or CC/TC ratio to elucidate the FECG phenomenon of cardiac size dependence. Furthermore, a comprehensive study that includes fetal ECG, Doppler and phase-rectified signal averaging measurements would potentially elucidate a more accurate cardiophysiologic model of IUGR.

Cardiac remodeling in IUGR is at least a three-dimensional process where anatomical changes cannot be considered alone but synchronously as a function of time and fetal size as well. In unaffected fetuses, the range of normal HR for instance varies with gestational age and can be observed as changing with developing maturity. This is because fetal heart structure and function continues to develop all throughout gestation and in the first year after delivery. This elasticity or penchant for adaptation is no longer a property of adult human hearts and is crucial in the configuration of diagnostic models in the fetus. For IUGR fetuses, setting cut off values for predictive markers will have to take multivariate algorithms into consideration. The cut off values for each ECG parameter in this study showed that fetal ECG is sensitive in detecting ongoing fetal cardiac remodeling but not a specific tool as it cannot be used to distinguish the cause of cardiac change, i.e., it cannot be used to distinguish between IUGR, congenital anomalies or even twin-to-twin transfusion syndrome. The potential for fetal ECG as a specific diagnostic tool lies in developing its vector analysis capability for determining specific heart structures involved (e.g., inferior heart wall, left ventricle or ventricular septum).

Detection of the subtle signs of cardiovascular remodeling is critical to monitoring hypertension and hypervolemia

#### REFERENCES


experienced by IUGR fetuses. The availability of noninvasive fetal ECG technology has ushered in the newest frontier of fetal cardiac studies to address this. It will increase our knowledge of this little-known area of electrocardiophysiology and ultimately, benefit our smallest of patients.

#### ETHICS STATEMENT

This research conformed with appropriate Good Clinical Practice (GCP) Guidelines. The protocol was approved by the respective research ethics boards of both the University of the Philippines-Manila [UPMREB (OBG) 2015-065-01] and the Harris Birthright Research Centre, Kings College, London (REC 02-03-033) followed by memorandums of agreement for collaborative research drawn between these centers prior to commencement of the study. All subjects gave written informed consent in accordance with the Declaration of Helsinki.

#### AUTHOR CONTRIBUTIONS

CV designed the study; all authors substantially contributed to the experiments; CV and KN handled all patient recruitment and care in London; CV and JS handled all patient recruitment and care in Manila. YK and KF extracted and analyzed the fetal ECG at Tohoku University; CV wrote the manuscript with contributions from all authors.

## FUNDING

This study was supported by grants from Grants-in-Aid for Scientific Research (09154228) and the Fetal Medicine Foundation (Charity no.: 1037116).

#### ACKNOWLEDGMENTS

The authors would like to acknowledge T. Ito and M. Endo for technical support; and M. Borja and J. Batara for statistical support.


recorded by the STAN and Nottingham systems. BJOG 101, 582–586. doi: 10.1111/j.1471-0528.1994.tb13647.x


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Velayo, Funamoto, Silao, Kimura and Nicolaides. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Hybrid EMD-Kurtosis Method for Estimating Fetal Heart Rate from Continuous Doppler Signals

Haitham M. Al-Angari <sup>1</sup> \*, Yoshitaka Kimura<sup>2</sup> , Leontios J. Hadjileontiadis 3, 4 and Ahsan H. Khandoker <sup>1</sup> \*

<sup>1</sup> Department of Biomedical Engineering, Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates, <sup>2</sup> The Graduate School of Medicine, Tohoku University, Sendai, Japan, <sup>3</sup> Department of Electrical and Computer Engineering, Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates, <sup>4</sup> Department of Electrical and Computer Engineering, Aristotle University of Thessaloniki, Thessaloniki, Greece

Monitoring of fetal heart rate (FHR) is an important measure of fetal wellbeing during the months of pregnancy. Previous works on estimating FHR variability from Doppler ultrasound (DUS) signal mainly through autocorrelation analysis showed low accuracy when compared with heart rate variability (HRV) computed from fetal electrocardiography (fECG). In this work, we proposed a method based on empirical mode decomposition (EMD) and the kurtosis statistics to estimate FHR and its variability from DUS. Comparison between estimated beat-to-beat intervals using the proposed method and the autocorrelation function (AF) with respect to RR intervals computed from fECG as the ground truth was done on DUS signals from 44 pregnant mothers in the early (20 cases) and late (24 cases) gestational weeks. The new EMD-kurtosis method showed significant lower error in estimating the number of beats in the early group (EMD-kurtosis: 2.2% vs. AF: 8.5%, p < 0.01, root mean squared error) and the late group (EMD-kurtosis: 2.9% vs. AF: 6.2%). The EMD-kurtosis method was also found to be better in estimating mean beat-to-beat with an average difference of 1.6 ms from true mean RR compared to 19.3 ms by using the AF method. However, the EMD-kurtosis performed worse than AF in estimating SNDD and RMSSD. The proposed EMD-kurtosis method is more robust than AF in low signal-to-noise ratio cases and can be used in a hybrid system to estimate beat-to-beat intervals from DUS. Further analysis to reduce the estimated beat-to-beat variability from the EMD-kurtosis method is needed.

haitham.alangari@kustar.ac.ae; haitham.angari@gmail.com

Specialty section:

Edited by: Raimond L. Winslow, Johns Hopkins University,

> United States Reviewed by: Amira J. Zaylaa,

United States \*Correspondence: Haitham M. Al-Angari

Ahsan H. Khandoker ahsan.khandoker@kustar.ac.ae

This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology

Lebanese University, Lebanon Martin Gerbert Frasch, University of Washington,

Received: 21 May 2017 Accepted: 15 August 2017 Published: 30 August 2017

#### Citation:

Al-Angari HM, Kimura Y, Hadjileontiadis LJ and Khandoker AH (2017) A Hybrid EMD-Kurtosis Method for Estimating Fetal Heart Rate from Continuous Doppler Signals. Front. Physiol. 8:641. doi: 10.3389/fphys.2017.00641 Keywords: fetal heart rate, fetal Doppler ultrasound, autocorrelation, kurtosis, empirical mode decomposition

# INTRODUCTION

Fetal Heart Rate (FHR) monitoring (1997) by Doppler based cardiotocography (CTG) in the third trimester is a commonly established method to identify fetal compromises. All pregnancies are usually checked by FHR monitor to identify any abnormality in FHR pattern (2001). The decrease of the FHR indicates an abnormal situation of the pregnancy particularly during uterine contraction (Peters et al., 2001). It has been shown however that modern fetal monitors using the Doppler US technique (DUS) do not provide reliable evaluation of FHR variability (Fuchs, 2014). This is due to the changing fetal Doppler signal over time as a result of location changes of fetal heart and the high signal to noise ratio (SNR) which make it very difficult to determine the beat-tobeat intervals (Shakespeare et al., 2001). Also, abdominal fetal electrocardiography (fECG) has been reported to provides more reliable description of the instantaneous FHR variability than DUS-based approaches due to the improvement in instrumentation, electrode technology and signal processing approaches related to detecting fECG from abdominal maternal ECG (Jezewski et al., 2017). Another limitation with the DUS techniques (in the currently available clinical systems) is that the estimated FHR typically resampled at 4 Hz which intrinsically dilutes the data stream of the short-term FHR fluctuations (Durosier et al., 2014; Li et al., 2015).

However, the important question still remains to be answered if it is possible to obtain beat-to-beat FHR reliably from fetal Doppler signals as compared to the same obtained by R–R intervals of fECG signals?

Various processing strategies have been tested on Doppler signals to generate FHR. For examples, band-pass filtering (Tuck, 1982; Spencer et al., 1987; Boos and Schraag, 1992) and the use of auto-correlation function (AF) (Tuck, 1982). However, it has been shown that the estimated heart rate variability (HRV) computed by the AF is inaccurate (Takeuchi and Hogaki, 1978; Divon et al., 1985). This is due to the averaging nature of the AF method where a single periodicity value is determined from all (neighboring) heart beats enclosed in the AF window (Lauersen et al., 1976; Cesarelli et al., 2009). Also, the AF uses the envelope of the DUS signal to determine the instantaneous periodicity which correspond to the cardiac cycle but not the consecutive heart beats. Recently an AF method with an adaptive window size has shown improved accuracy in estimating the beat-to-beat interval from DUS (Jezewski et al., 2011). Using that method, the cardiac cycles are represented by sets of periodicity measurements (computed from the autocorrelation of multiple shifts of a certain window size). Then, these periodicity measurements are segmented using a segmentation algorithm where the location of the successive segment is continuously adjusted and based on the measurements contained in each segment an estimated of the new cardiac cycle duration is computed.

The empirical mode decomposition (EMD) method decomposes a signal into components with well-defined instantaneous frequency. These components are called Intrinsic Mode Functions (IMFs). Each IMF has a unique local frequency and different IMFs do not exhibit the same frequency at the same time (Huang et al., 1998). The ensembled EEMD evolved from EMD to treats the problem of mode mixing where different modes of oscillation may appear in one IMF or one mode can spread across different IMFs (Wu and Huang, 2009). This make the EEMD a true filter for any data. Applying the kurtosis on the (filtered) IMFs is done to determine which parts of these IMFs is important (high kurtosis value) and unimportant (low kurtosis value). Kurtosis has been used previously in detecting waveform changes (Saragiotis et al., 2004; Rekanos and Hadjileontiadis, 2006). In the DUS case, the locations of the signal peaks are expected to be within these high kurtosis parts.

Therefore, in this study, we propose a method combining EMD and kurtosis to estimate FHR and fetal HRV from DUS signal. We compare FHR (mean, standard deviation, and root mean squared successive difference) values estimated from DUS using both the AF method and our proposed method with a true RR interval determined from fECG of 44 healthy fetuses from early and late gestational age groups.

# MATERIALS AND METHODS

#### Subjects

Abdominal ECG and Doppler ultrasound signals (DUS) were collected from 51 pregnant women at Tohoku University Hospital in Japan. The pregnant women were lying on their backs while the abdominal ECG signals were collected using 12 electrodes: ten on the mother's abdomen, one reference electrode on the back and one electrode at the right thoracic position. Signals were recorded during daytime (between 10 a.m. and 4 p.m.) over 3 years (2009–2011). The same experimental set up was applied to all the pregnant mothers who participated in this study. The continuous DUS data were obtained using ultrasonic transducer 5,700 (fetal monitor 116, Corometrics Medical Systems, Inc.) with 1.15 MHz signals. The recordings were of 1-min length and were sampled at 1 kHz with 16-bit resolution. The continuous DUS was recorded in a laptop by a data acquisition system which synchronizes with Abdominal ECG machine. fECG signals were separated by another custom-made software. The study protocol was approved by Tohoku University Institutional Review Board (IRB: 2015-2-80-1) and written informed consent was obtained from all subjects. The inclusion criteria for the study were: (1) Signed on written consent form, (2) Maternal age of 20 years or older, and (3) gestational age in the range of 24– 42 weeks. The exclusion criteria were: (1) Diagnosed with multiple pregnancy, abnormal pregnancy, pregnancy with an obstetric complication (e.g., gestational diabetes, gestational hypertension, uterine fibroids, and cervical cancer) and (2) Scheduled for Caesarean section. The raw ECG and DUS signals were visually checked and noisy records (visually no peaks were seen) were removed from the dataset (two records). The remaining records were divided into two age groups: early gestational group (≤32 weeks; 20 cases) and late gestational group (≥35 weeks; 24 cases). Records with gestational age inbetween the two groups were removed from the analysis (five records).

# Empirical Mode Decomposition (EMD)

In this work, the EMD was used to analyzed the oscillatory behavior of DUS. The EMD method deals with non-stationary and non-linear data (Huang et al., 1998). The EMD method consider that a signal consists of different simple intrinsic modes of oscillations. These simple oscillations are represented by the IMFs which satisfy two conditions:


A signal x(t) is represented by the EMD method as:

$$\mathbf{x}(t) = \sum\_{i=1}^{N} \mathbf{C}\_i(t) + r\_N(t), \tag{1}$$

where ci(t) is the i-th IMF and rN(t)the final residue.

#### Estimation of R Peak Location

To separate fECG from the composite abdominal signal, a combination of maternal ECG cancelation and blind source separation with the reference signal (BSSR) was used (Sato et al., 2007). In brief, electrical activities of the heart can be modeled as a vector in the direction of excitation called the heart vector (Symonds et al., 2001). The maternal ECG component was excluded by subtracting the linear combination of mutually orthogonal projections of the heart vector. After that, BSSR which is a kind of neural network method, was used to extract fECG from complex mixture using DUS signal as a reference. RR interval (intervals between successive R waves of the fECG were then computed using the algorithm developed by Pan and Tompkins (1985). The computed RR intervals from the fECG (true RR) were considered, the ground truth to be compared with the estimated beat-to-beat intervals from the DUS using the combined EMD-kurtosis and the AF methods.

To estimate beat location from DUS using the combined EMD-kurtosis method, the background noise was first removed from the DUS signal using wavelet thresholding with the first 15 levels of Haar wavelet. The signal was then decomposed into its IMFs using the EEMD method.

The kurtosis defined as:

$$\hat{\wp}\_4 = (N-1) \frac{\sum\_{n=1}^N \mathfrak{x}^4(n)}{\left(\sum\_{n=1}^N \mathfrak{x}^2(n)\right)^2} \tag{2}$$

was computed for each IMF signal x(n) where N is the number of the sample in the signal. The kurtosis of each IMF was tested against the Chebyshev inequality to determine whether it was good or noisy for the detection of the DUS peak location. More details on the method can be found at (Papadaniil and Hadjileontiadis, 2014).

Then, for all selected IMF a sliding window with a width of 50–600 ms (with an increment of 50 ms) and a shift of 1 ms was used to compute the kurtosis at each shift point. A matrix of IxW kurtosis vectors was constructed where I indicate the number of selected IMF and W is the number of used windows. Again, each vector was tested against the Chebyshev inequality to be selected in the final subgroup used to estimate beat location. Two other measures were introduced to fine select the best kurtosis vectors:

(1) The percentage of the mismatch error between true R peaks count (determined from fECG) and estimated beat count (determined by the peaks of the selected kurtosis vectors) as:

$$mismatch\,error = \frac{true\,\,R\,count - estimated\,\,heat\,\,count}{true\,\,R\,count} \times 100\,\,\,\,\,\tag{3}$$

(2) The variability of the estimated beat location (varB\_loc) was determined by the standard deviation of the absolute difference between real R peak location and the estimated beat location:

A selected kurtosis vector would minimize both the mismatch error and varB\_loc for better estimation of the beat location. Finally, the selected vectors were summed up and a peak search with minimum peak distance of 300 ms was performed.

The estimation of beat-to-beat intervals using AF method was done similar to (Jezewski et al., 2011). Two measures were used to compare results between the two methods. These were:


FIGURE 1 | An example of fetal Doppler signal and the derived kurtosis signal: (A) 5-s segment of fetal ECG (fECG), (B) corresponding Doppler ultra sound signal (DUS), (C) DUS after noise removal using the first 15 levels of the Haar wavelet, (D) constructed kurtosis signal from the selected IMFs of DUS (green stars are location of the kurtosis peaks and red stars are true location of the R peaks determined from fECG). Signals were normalized to their maximum values.

$$\begin{aligned} \text{Successive} \quad &\text{beat } error \text{ (i)} = \\ &\left| \frac{true \text{ RR } (i) - estimated \text{ beat} - to - beta (i)}{true \text{ RR} (i)} \right| \times 100 \end{aligned} \tag{4}$$

for i = 1: Nmin where Nmin is minimum length of the two true RR and estimated beat-to-beat vectors.

Finally, a comparison between the true RR and estimated beat-to-beat (using AF and EMD-kurtosis methods) in the mean, standard deviation (SDNN) and root mean square of successive difference (RMSSD) was performed using Kruskal-Wallis test. Significant differences were reported if p < 0.05.

#### RESULTS

The mean ± Standard Deviation (SD) of the early group gestational age was 27 ± 4.3 weeks and for the late group was 38 ± 1.7 weeks. The RR interval duration of the early group was 406.4 ± 27 ms and for the late group was 422.3 ± 29 ms. An example of an extracted fECG with the corresponding DUS and the estimated beat location using the EMD-kurtosis method is shown in **Figure 1**. After noise removal from the raw DUS (**Figures 1B,C**), the IMFs were generated and the kurtosis of the IMFs with different window size was computed. The kurtosis signals with optimum window sizes and EMD levels

FIGURE 2 | Absolute value of mismatch error and varB\_loc measure used to optimize selection of window size to compute the EMD-kurtosis signal. For the early group the minimum absolute mismatch error and varB\_loc were achieved at window size between 300 and 350 ms (A,C, respectively) and for the late group the minimum absolute mismatch error and varR\_loc were achieved at window size between 350 and 400 ms (B,D, respectively). The arrow point at window sizes that minimize the measurements.

TABLE 1 | Comparison between the EMD-kurtosis and the AF methods in mismatch error and mean successive beat error.


\*,† AF is significantly higher than EMD-kurtosis with p < 0.05, p < 0.01, respectively.

‡ values represented in mean root squared. were summed up (**Figure 1D**) and the peaks of the resultant signal were detected. The differences between the successive peak locations represented the estimated beat-to-beat intervals. The optimum EMD levels and window sizes were selected based on the criteria shown in **Figure 2**.

**Figure 2** shows a surf plot of the average mismatch error and varB\_loc for the early and late age groups. For the early group, the minimum mismatch error and varB\_loc were achieved at window size around 300–350 ms and the first three EMD levels while for the late group the minimum was achieved around window size of 350–400 ms with same EMD levels. This smaller grid was used to compute the final (summed) kurtosis signal.

**Table 1** shows a quantified comparison between the EMDkurtosis and AF methods. Both the mismatch and mean successive beat errors were lower using EMD-kurtosis for both early and late groups. This can also be observed from **Figure 3**. The mismatch error was significantly lower using EMD-kurtosis for the early group (p < 0.01) and the mean successive beat error was significantly lower using EMD-kurtosis in the late group (p < 0.05). No significant was found when comparing early and late groups using the same method. **Figure 3** shows the % mismatch and mean absolute beat-to-beat error for all cases in the early and late groups. Generally speaking, these measures were lower than the 10% error limit for the EMD-kurtosis method while for the AF method some case exceeded that limit. % mismatch error was clearly higher using AF in the early group (**Figures 3A,B**). Four cases in the late group were having higher absolute beatto-beat error compared to the other cases in the group using

TABLE 2 | Comparison between true RR and estimate beat-to-beat intervals using the EMD-kurtosis and the AF methods.

AF method (B,D). Negative sign of mismatch error indicates more estimated beats than the true value.


† Significantly different from true RR (p < 0.01). Kruskal-Wallis test was done between true RR parameters and estimated parameters from each of the EMD-kurtosis and the AF methods separately.

AF which was not observed using the EMD-kurtosis method (**Figures 3C,D**).

**Table 2** shows a comparison of the mean, SDNN and RMSSD between the true RR and estimated beat-to-beat. The EMDkurtosis method was better in estimating mean beat-to-beat while the estimated SNDD and RMSSD were significantly higher than the true RR values using both methods (the EMD-kurtosis showed higher values).

**Figure 4** helps interpreting these results with examples of true RR and estimated beat-to-beat intervals for an early and a late case. When true RR variability is small the AF was better in estimating beat-to-beat intervals (**Figures 4A,B**) while the EMD-kurtosis method worked better with high RR variability.

Bland-Altman plots of the estimated mean beat-to-beat, SDNN and RMSSD using the two methods compared to the true RR interval computed from fECG are shown in **Figure 5**. The EMD-kurtosis method has smaller mean difference between the true RR and the estimated beat-to-beat intervals mean where most of the cases were within the mean ± 2 SD range while the AF has smaller mean difference between true and estimated SDNN and RMSSD.

#### DISCUSSION

In this work, we have applied a method based on computing the kurtosis function on the IMFs extracted from the DUS signal to estimate cardiac beat-to-beat intervals. The proper choice of the window size used to compute the kurtosis was important in reducing the error in estimating both the number beats and the beat-to-beat variability. A window size slightly smaller than the duration of the mean RR interval provided the optimum results. This is obvious since a window of this size results in almost one distinguished kurtosis peak in a single cardiac cycle. Smaller or bigger windows increase the changes of missing true peaks in the situation of low SNR (the early group) and of having multiple peaks due to distinct walls and valves movements (the late group).

Due to developing fetal heart wall and valves in the early gestational week, the reflected Doppler signal is relatively weaker which reduces the SNR. This could be the reason of having higher error for AF method in the early gestational group (even when compared with the results of the same method for the late group). On the other hand, the EMD-kurtosis method was more robust against low SNR. It reduced the mismatch error four times for the

method showed better estimation (A,B) while the EMD-kurtosis method performed better for signals with higher RR variability (C,D).

early group and two times for the late group as compared to the AF method.

The mean successive beat error for the EMD-kurtosis method although was lower than AF but was still considered high (above 5%, **Table 1**) for good evaluation of HRV. Various reasons could explain this high error. Firstly, the DUS represents the mechanical movement of heart walls and valves while the R wave represents the electrical depolarization of the ventricle which results in an expected delay between the two signals. This delay could vary from beat to beat.

Also, heart valves movements which considered of high frequency cannot be completely separated from wall movements which have lower frequency which results in an increased variability in estimating HRV from the DUS signal. Other factors including movement of mother or baby and changing the orientation of the Doppler probe to fetal heart could also cause variability in estimating HRV.

Our proposed EMD-kurtosis method showed also improvement in estimating mean heart rate compared to AF in both early and late groups (**Table 2**). The beat-to-beat variability assessed by SDNN and RMSSD was significantly higher than the true RR intervals for both methods although it was higher for the EMD-kurtosis method. This could be due to the adaptive AF window that changes its size slowly which on one hand helps in reducing RR variability (**Figure 4B**) but on the other hand responds slowly to rapid changes in the true signal (**Figure 4D**) which results in low estimated mean.

A limitation of our study is that the collected data was of 1-min length which is too short to control the fetal states. However, short term recording is typically used for clinical investigation in pregnancy clinic. Changing the segment length from 5 to 2 min has shown changes in the RMSSD values for FHR estimated by fetal magnetocardiogram (fMCG) which is known to have higher temporal resolution than the DUS signal (Moraes et al., 2012). Longer signal durations (10–30 min) will be needed to verify these findings. The EMD-kurtosis algorithm is a time-domain method that can easily be implemented in a microprocessor in the same DUS machine or in a separate device for practical clinical use. One consideration will be the time needed to extract the IMFs especially with long signal duration. This work is considered a step toward good estimation of FHR and its variability from DUS signal. Further research

#### REFERENCES


combining the EMD-kurtosis with other non-linear signal processing methods are needed to reduce variability in beat-tobeat intervals estimation and improve the overall accuracy.

#### AUTHOR CONTRIBUTIONS

HA, YK, LH, and AK designed the study approach. YK and AK designed the experiments. YK supervised data collection. HA was responsible for analysis and running statistics. HA, LH, and AK interpreted the results and wrote the manuscript.

### FUNDING

This study was supported by research incentive funds from Khalifa University Internal Research Fund Level 2.

and kurtosis features. IEEE J. Biomed. Health Inform. 18, 1138–1152. doi: 10.1109/JBHI.2013.2294399


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Al-Angari, Kimura, Hadjileontiadis and Khandoker. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Is Abdominal Fetal Electrocardiography an Alternative to Doppler Ultrasound for FHR Variability Evaluation?

Janusz Jezewski <sup>1</sup> , Janusz Wrobel <sup>1</sup> \*, Adam Matonia<sup>1</sup> , Krzysztof Horoba<sup>1</sup> , Radek Martinek <sup>2</sup> , Tomasz Kupka<sup>1</sup> and Michal Jezewski <sup>3</sup>

*1 Institute of Medical Technology and Equipment ITAM, Zabrze, Poland, <sup>2</sup> Department of Cybernetics and Biomedical Engineering, VSB-Technical University of Ostrava, Ostrava, Czechia, <sup>3</sup> Institute of Electronics, Silesian University of Technology, Gliwice, Poland*

#### Edited by:

*Ahsan H. Khandoker, Khalifa University, UAE*

#### Reviewed by:

*Haitham Alangari, Khalifa University, UAE Faezeh Marzbanrad, Monash University, Australia Yoshitaka Kimura, Tohoku University, Japan*

> \*Correspondence: *Janusz Wrobel januszw@itam.zabrze.pl*

#### Specialty section:

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology*

> Received: *20 March 2017* Accepted: *27 April 2017* Published: *16 May 2017*

#### Citation:

*Jezewski J, Wrobel J, Matonia A, Horoba K, Martinek R, Kupka T and Jezewski M (2017) Is Abdominal Fetal Electrocardiography an Alternative to Doppler Ultrasound for FHR Variability Evaluation? Front. Physiol. 8:305. doi: 10.3389/fphys.2017.00305* Great expectations are connected with application of indirect fetal electrocardiography (FECG), especially for home telemonitoring of pregnancy. Evaluation of fetal heart rate (FHR) variability, when determined from FECG, uses the same criteria as for FHR signal acquired classically—through ultrasound Doppler method (US). Therefore, the equivalence of those two methods has to be confirmed, both in terms of recognizing classical FHR patterns: baseline, accelerations/decelerations (A/D), long-term variability (LTV), as well as evaluating the FHR variability with beat-to-beat accuracy—short-term variability (STV). The research material consisted of recordings collected from 60 patients in physiological and complicated pregnancy. The FHR signals of at least 30 min duration were acquired dually, using two systems for fetal and maternal monitoring, based on US and FECG methods. Recordings were retrospectively divided into normal (41) and abnormal (19) fetal outcome. The complex process of data synchronization and validation was performed. Obtained low level of the signal loss (4.5% for US and 1.8% for FECG method) enabled to perform both direct comparison of FHR signals, as well as indirect one—by using clinically relevant parameters. Direct comparison showed that there is no measurement bias between the acquisition methods, whereas the mean absolute difference, important for both visual and computer-aided signal analysis, was equal to 1.2 bpm. Such low differences do not affect the visual assessment of the FHR signal. However, in the indirect comparison the inconsistencies of several percent were noted. This mainly affects the acceleration (7.8%) and particularly deceleration (54%) patterns. In the signals acquired using the electrocardiography the obtained STV and LTV indices have shown significant overestimation by 10 and 50% respectively. It also turned out, that ability of clinical parameters to distinguish between normal and abnormal groups do not depend on the acquisition method. The obtained results prove that the abdominal FECG, considered as an alternative to the ultrasound approach, does not change the interpretation of the FHR signal, which was confirmed during both visual assessment and automated analysis.

Keywords: Doppler ultrasound, fetal electrocardiogram, fetal heart rate analysis, fetal state assessment, fetal outcome

# INTRODUCTION

Fetal heart activity is a primary source of information which enables assessment of the fetal state during pregnancy and at labor. This information is obtained mainly through analysis of the fetal heart rate (FHR) signal being formed from the instantaneous values calculated according to the formula: FHR [bpm] = 60000/T [ms]. The FHR values are expressed in beats per minute, and T is the time interval between two consecutive fetal heart beats that comprises one complete cardiac cycle. Together with additional signals describing the uterine contractile activity and fetal movement profile, the FHR signal constitutes the cardiotocographic record. Acquisition of these additional signals is quite simple, but measurement of the fetal heart rate has been always a challenge. Already in 1960s the fetal electrocardiogram was recorded for the first time by means of electrode attached to fetal head. The quality of such recorded direct fetal electrocardiogram (FECG) is usually very good, and thus it enables, using a quite simple processing method, determination of the beat-to-beat intervals with very high accuracy. However, the invasive approach and application limited to the labor only caused that direct method did not found wide application in clinical practice.

As a result of further research and development the noninvasive ultrasound (US) method has become a standard approach since early 1970s, as it can be used both during pregnancy and labor. At present, all bedside fetal monitors intended to use in clinical conditions are based on the pulsed Doppler ultrasound technique, with measurement transducer attached to maternal abdomen. Principle of operation relies on internal processing of the envelope of the US beam reflected from moving parts of fetal hearts—valves or walls, to find the episodes corresponding to consecutive heart beats. However, a complex structure and varying content of the US signal, usually caused by relocation of the fetal heart in relation to a transducer during monitoring session, make a determination of the beat-to-beat interval very difficult (Khandoker et al., 2009; Marzbanrad et al., 2014). Therefore, a correlation techniques, considering full shape of the analyzed signal, have been applied. The cross-correlation technique with changeable template appeared to be too sensitive to US signal changes, which resulted in considerable signals loss. Thus, an autocorrelation function with adaptive window selection has been applied in next generation fetal monitors. However, the autocorrelation function does not detect the consecutive heart beats but only determines the instantaneous periodicity of the US signal envelope which corresponds to cardiac cycle being measured. This leads to effect of averaging of neighboring cardiac cycles and thus decreasing of FHR determination accuracy in relation to fetal electrocardiography (Lee et al., 2009; Voicu et al., 2010). The obtained FHR signal is provided by the bedside monitor as the trace in a printout with established time scale of 1, 2, or 3 cm/min. As long as the FHR trace has been analyzed visually, the lower accuracy did not affect significantly the fetal state assessment. More important was to ensure the trace continuity which allowed clinicians to observe a general tendency of the fetal heart rate changes, and to recognize the features representing longitudinal FHR patterns relating to the fetal state, like acceleration or deceleration. It was found that the evaluation of fetal state, when based on visual interpretation, has been mainly affected by low inter- and intraobserver agreement (Jezewski et al., 2002; Romano et al., 2016a). That was a result of both complexity of the FHR signal and the fact that important part of information relating to instantaneous changes of FHR values has been hidden from a naked eye. These changes are considered to be very important FHR characteristics, reflecting appropriate neurological modulation of the FHR.

Thus, further development stage of fetal monitoring was aimed at automated analysis of the FHR signal and its implementation as built-in procedure of bedside monitors as well as in computer-aided fetal monitoring system. Some other requirements important for monitoring the pregnant women in hospital, like surveillance of many patients, detecting and alerting of symptoms of fetal distress, or electronic archive with the signals and perinatal data, have made the computer-aided system with online automated analysis the standard in modern obstetrics (Wrobel et al., 2013, 2015a). Automated analysis comprises detection and description of the above mentioned FHR features, like acceleration and others, as well as determination of the instantaneous FHR changes by providing a set of indices to evaluate the long-term and short-term (beat-to-beat level) variability of the fetal heart rate.

Automated online analysis provides a quantitative description of the FHR, but the final interpretation of the record is still done by a clinician. There are a number of papers relating to automated classification of the FHR recordings by using different methods of computational intelligence like neural network, support vector machines or epsilon-insensitive learning (Czabanski et al., 2008, 2013). However, taking into account that the input data set comprised the automatically determined features of the FHR signals, collected from the clinical databases (Chudacek et al., 2014), the obtained classification results should be faced to the limitations of the ultrasound approach as it is discussed below (Voicu et al., 2014; Wrobel et al., 2015c).

The variability indices were originally defined using the beatto-beat intervals determined from the direct electrocardiogram. Their straight application to the FHR signals being provided by fetal monitors raised a question how a limited accuracy of the ultrasound approach affects determination of the cardiac intervals, and thus the variability indices values. Several research studies were aimed at evaluation of the reliability of the ultrasound method in reference to the direct electrocardiography (Ibrahimy et al., 2003; Reinhard et al., 2010, 2012; Cohen et al., 2012; Kimura et al., 2012; Desai et al., 2013; Kording et al., 2015). However, they were aimed at comparing the signal loss episodes or directly the FHR values. In our study we showed that the error of cardiac cycle determination (instantaneous FHR value) has not been correlated with the FHR variability indices error. It means that the measurement accuracy resulting from the fetal monitor specification cannot be directly related to the results of the computer-aided analysis of the FHR variability. In general, we concluded that modern fetal monitors using the Doppler US technique are not able to provide the signal with the accuracy required for reliable quantitative evaluation of instantaneous FHR variability, particularly the short-term variability, based on the indices calculated automatically (Voicu et al., 2014). Fortunately, the values of indices determined in that way are underestimated, which prevents the fetal distress signs from being undetected.

Several attempts were carried out to improve the reliability of the ultrasound method by using advanced signal processing of the Doppler envelope (Jezewski et al., 2011), but none of them have been applied in the bedside monitors yet. In Wrobel et al. (2008) the method has been proposed to improve the reliability of the FHR variability indices, which relies on the errors recognized in the ultrasound measurement channel.

The fact, that the FHR signal obtained from the ultrasound approach has been recognized as not good enough to fully exploit the potential of automated analysis offered by computeraided fetal monitoring system, brought back an interest of the fetal electrocardiography (Fuchs, 2014). However, taking into account the need to monitor a whole pregnancy period, only a noninvasive approach could be considered, which relies on indirect recording the FECG from electrodes located on maternal abdominal wall (Ungureanu et al., 2009; Vullings et al., 2010; Kimura et al., 2012; Khalaf et al., 2013; Behar et al., 2014; Agostinelli et al., 2015b).

Another important issue for development of effective abdominal electrocardiography refers to a growing interest in high-risk pregnancy telemonitoring at home (Wrobel et al., 2015b). When using the ultrasound-based fetal monitor the transducer has to be carefully placed to ensure the ultrasound beam is focused on the fetal heart. What's more, during monitoring session the transducer may require repositioning due to a change of fetus position. Otherwise, the signal loss occurs which may cause, in case when a woman performs the monitoring session alone, her unfounded fear and unpredictable reaction. When the abdominal electrodes are fixed on the abdomen, the patient can easily verify the signal loss, which in that case occurs only when one of the electrodes peels off (Karvounis et al., 2007; Kolomeyets and Roshchevskaya, 2013; Agostinelli et al., 2015b).

Improvement of the measurement instrumentation, electrode technology and the signal processing methods that have been noticed during recent years, enabled to cope with the problems connected with development of the abdominal fetal electrocardiography (Kotas, 2008; Vullings et al., 2010). The signal acquired from fetus head is in fact "pure" fetal electrocardiogram, whereas the abdominal signal includes also the maternal electrocardiogram (MECG) and some noise coming mainly from muscle activity (Taralunga et al., 2009, 2014; Martinek et al., 2016). Thus, the crucial step in extraction of the FECG from the abdominal signal is a suppression of maternal electrocardiogram while preserving the fetal QRS complexes (Melillo et al., 2014; Agostinelli et al., 2015a). The energy of MECG is many times higher than the energy of FECG, and what's more the frequency band of both these components partly overlaps which makes simple filtering useless (Karvounis et al., 2007). A number of different approaches to MECG suppression and detection of fetal QRS complexes were presented in literature (Ungureanu et al., 2007; Liu and Luan, 2015; Poian Da et al., 2016). The system for acquisition of abdominal signals and original method for FECG extraction were proposed by the authors, and the indirect fetal electrocardiography was evaluated in relation to the gold standard—direct FECG approach (Jezewski et al., 2012). Referring to the results obtained in our previous study, concerning a comparison of ultrasound approach with the direct FECG (Jezewski et al., 2006), we concluded that the abdominal fetal electrocardiography provides accuracy not worse than the ultrasound method does. However, in all studies where the US method or abdominal FECG was compared with direct FECG, the results were obtained only for the signals being acquired during labor. Considering that fetal development, taking place during a whole pregnancy period, affects the characteristics of the FHR signals, we decided to carry out the comparison of the abdominal FECG and ultrasound method based on the signals collected during pregnancy. It is obvious that such approach excludes the direct electrocardiography from the study, and causes some problems for comparison methodology due to a lack of reference data (Sato et al., 2011; Cohen et al., 2012; Kimura et al., 2012). Since both types of the signals were acquired by means of two sets of instrumentation, another important problem has been recognized—the FHR signals synchronization, i.e., finding the corresponding cardiac cycles. It should be noticed that in case of ultrasound-based monitor, the FHR signal is provided through its output only as the measurement values of instantaneous heart rate evenly spaced with 250 ms. On the other hand, the system for noninvasive FECG is able to provide, along with the evenly spaced signal, the time event series with durations of consecutive cardiac cycles.

In this work the methodology is proposed to compare two different methods for fetal heart rate monitoring. Its originality relates to the fact that comparison has been carried out not only in relation to the corresponding cardiac cycle values, but also to the clinically important indices describing the instantaneous FHR variability.

### METHODS

The research material comprised the FHR signals acquired simultaneously using the Doppler ultrasound as well as the electrocardiographic methods in a group of 70 pregnant women. From a number of monitoring sessions performed for each patient, we selected only one recording acquired around 1 week before delivery, with a length of at least 30 min (Georgieva et al., 2014). All the recordings are accompanied by information on fetal outcome: gestational age at birth, blood gas parameters pH and BE, percentile of fetal birth weight, Apgar score, information about a possible stay in the NICU. The patients were monitored by simultaneously using two popular in maternity wards, systems for fetal and maternal monitoring: MONAKO and KOMPOREL. Unfortunately, these systems were unable to synchronize recorded signals during the monitoring session. The time shift between signals beginnings in each session could reach up to a few minutes, whereas in case of comparative studies the precise synchronization is required (even on the level of individual heartbeats). Hence, the problem of signals synchronization has been considered as a significant challenge.

As the result of each simultaneous monitoring session, two files were obtained of the native format, where the FHR signal is represented by the values measured evenly with 250 ms period. Files from the MONAKO System comprise the FHR\_U signal captured from the output of fetal monitor equipped with the ultrasound transducer (Hewlett-Packard M1351). The files from the KOMPOREL System provide the FHR\_E signal being determined on a basis of fetal electrocardiogram recorded from the abdominal wall of the mother. The fetal electrocardiogram is recorded by using four electrodes placed on the maternal abdomen. The crucial step in extraction of the FECG from the abdominal signal is a suppression of maternal electrocardiogram while preserving the fetal QRS complexes (Castilloa et al., 2013; Martinek et al., 2015). The proposed method for the MECG suppression is based on subtracting the pattern of maternal P-QRS-T complexes and spatial filtering. It ensures correct determination of the fiducial points as well as the factors scaling the pattern. The algorithm for detection of the fetal QRS complexes is based on a matched filtering approach in order to reduce the sensitivity to interferences. Additionally, the detection is carried out with a set of decision rules to predict the duration of the next beat-to-beat cardiac cycle (Matonia et al., 2006). As a result of the FECG analysis, the time event series is obtained as the time markers when the successive fetal heartbeats were detected—the QRS complexes, which is then used to determine the FHR\_E signal as the values with 250 ms period (Guerrero Martinez et al., 2006; Almeida et al., 2013, 2014).

#### Signals Synchronization

The procedure for synchronization of each pair of the FHR\_U and FHR\_E signals consisted of two stages. In the first stage an initial visual adjustment was supported by a dedicated program for visualization of the signals. This program as well as all the others, created for the purpose of this work, was developed in LabView environment (NationaI Instruments) (Desai et al., 2013). After coarse synchronization of the signals the common part of FHR\_U and FHR\_E signals was separated. It relied on moving the beginning and end of one signal, to indicate the fragment of interest according to the other signal. At this stage the signals quality had to be good enough to allow recognizing the characteristic features common for both signals, and constituting the so called centering points, and the common part of the signals had to have non-zero length. These conditions were not met in case of three recordings, hence in further processing only the set of 67 patients were included.

In the second stage the signal validation was conducted, as well as precise synchronization of both signals, at the level of individual FHR values provided every 250 ms. The developed software enabled semi-automated synchronization. The program automatically found the time shift between the signals to ensure the minimum differences between the corresponding values, which mostly led to proper synchronization. After that, the visual verification was carried out with a possibility of additional time shift correction, followed by the final acceptance (**Figure 1**).

The software for determination of the optimal time shift between the analyzed signals was using the synchronization function based on the mean absolute differences (MAD), determined for the corresponding (applying the time shift) FHR\_U and FHR\_E values. To improve the performance of synchronization function it was necessary to further reduce the influence of random interferences appearing in the FHR signals, as well as sudden value changes resulting from the measurement errors or potential acceleration and deceleration episodes. Hence, the segments with sudden changes in the FHR signal were excluded from the function determination, if the absolute difference between a given value and the preceding one was higher than 10 bpm (Spilka et al., 2012). If a given FHR value was rejected, the next one was compared to the mean calculated from the previous values (including the rejected) in the 240 values window. The FHR values were also rejected from the signal, which were suspected to represent the maternal heart rate—the details of the algorithm are presented in Wrobel et al. (2015b). This type of erroneous measurements occurs in the FHR\_U signal, as it is typical for the ultrasound method and quite frequent in the US-based fetal monitors. The above-mentioned preprocessing is only intended to synchronize the signals and do not change their information content.

The optimum time shift, corresponding to the fully synchronized signals, was obtained for the function minimum, when applying additional shift in the range from −25 to +25 FHR values (measured every 250 ms), in respect to the signals synchronized after the first stage. If at the beginning or end of a given signal any interference associated with the start or end of the monitoring session occurred, they were also removed in the trimming process. Trimming to the full minutes in turn, results from the fact that the analysis of the instantaneous FHR variability is always carried out within a 1-min signal segments. Finally, as a result of the synchronization procedure some signal pairs could be shortened by as much as 4 min. After the second stage of synchronization the common part of the analyzed signals is trimmed to the largest whole number of the minutes (the number of FHR values was a multiple of 240).

We assumed that the minimum length of synchronized signals (constituting the pair) subjected to further analysis should be 10 min. According to that criterion only one recording was rejected, and 66 recordings were left.

#### Signal Loss Analysis

The next processing step consisted of verifying the signal pairs in terms of their quality, measured both by a size and nature of the signal loss episodes. Episodes of signal loss are preliminarily detected by the monitoring systems used, and represented by zero values in the FHR signal. For the purpose of this work also the potential erroneous FHR values, as not meeting the adopted criteria, were marked as signal loss episodes, using the dedicated developed program. Signal segments were considered to be signal loss if they did not meet the van Geijn modified criterion (van Geijn, 1980), proposed in Jezewski et al. (2012). That was applied to these FHR values, which were considered as errors by a procedure for sudden changes removal used in the second stage of synchronization. In this case, the established thresholds are: the absolute difference between a given FHR value and the preceding one greater than 20 bpm, and window width equals to 100 values. The segments, suspected to contain

the maternal heart rate signal, were indicated as the signal loss episodes—similarly as in the second stage of synchronization. Finally, as a result of the signal loss analysis, the FHR\_U and FHR\_E signals were obtained with additional information about detected gaps (FHR values equal to zero). The signal loss level is defined as a percentage of the duration of signal loss episodes (the number of FHR values equal to zero) in relation to the total duration of the signal (all FHR values). Taking into account the maximum level of signal loss of 30% in either FHR\_U or FHR\_E signal, five recordings were removed. Additionally, in terms of uniformity of signal loss distribution in time, four questionable recordings with signal loss between 20 and 30% were visually assessed. Only one recording was excluded due to the accumulation of the signal loss (equal to 23%) in the middle part of the FHR\_U signal. The final research material consisted of signal pairs from 60 monitoring session. The total length of recordings was equal to 1995 min. The length of individual recordings varied from 11 to 64 min, with an average of 33.3 min.

#### Direct Signals Comparison

For the final set of recordings, consisting of the FHR\_U and FHR\_E signal pairs, some descriptive statistics of signal comparison were calculated, both on a global basis as well as at the level of particular FHR values. These statistics include (calculated for individual recordings): the recording duration, signal loss level, mean value of the differences between the corresponding instantaneous FHR values (MD), standard deviation (SD), mean absolute difference (MAD), as well as the summary statistics for the entire research material. The difference between the pairs of corresponding instantaneous values of FHR\_U and FHR\_E was expressed dually: as the heart beats per minute (bpm) as well as in milliseconds. The second representation is obtained by conversion of the FHR values into intervals between successive heart beats, according to the hyperbolic transformation with the 60,000 factor.

As equally important it is assumed the comparison of the signals in a format commonly used in automated analysis of the low variability FHR signal components (Jezewski et al., 2002). In this format the successive values of the signal are determined by averaging 10 consecutive original FHR values (time series measured every 250 ms). Important part in the averaging process is the removal of the signal loss episodes, marked as zero values. If, while averaging 10 original FHR values, more than four values are marked as signal loss, the resulting value of 2.5 s period is also considered as signal loss and is assigned with zero value.

#### Indirect Signal Comparison via Clinical Parameters

As shown in the literature (Jezewski et al., 2016), the differences from the direct comparison of signals are often not correlated with differences in values of clinically important parameters of quantitative description of FHR signal, determined by the fetal monitoring systems (Georgieva et al., 2012). These parameters are used by clinicians, interpreting the FHR signals in order to assess the fetal state. Therefore, it was considered as important to identify the impact of the FHR signal acquisition method on the clinically significant parameters. For that purpose, the synchronized FHR\_U and FHR\_E signals were saved into the native format files (measured with 250 ms) and reloaded to the archive of MONAKO System. Thanks to an option of reanalysis of the archival records, for each of 60 recordings (120 FHR signals), the quantitative parameters describing the variability patterns detected in the FHR signal were determined automatically. As the result, the lists of parameters were obtained for the FHR\_U and FHR\_E signals. They included the parameters describing some patterns of FHR variability in the time domain: mean value of FHR (M\_FHR), mean value of the FHR baseline (M\_BL), number of detected acceleration (ACC) and deceleration (DEC) episodes (Georgieva et al., 2012; Wrobel et al., 2013). Additionally, an assessment of the instantaneous FHR variability was provided as: duration and value of the high (HE\_D and HE\_V) and low (LE\_D and LE\_V) variability episodes, average value of the longterm variability (LTV) and short-term variability (STV) indices, as well as the FHR oscillations (OSC) together with the percentage of different oscillation types (OSC\_I÷OSC\_IV) (Jezewski et al., 2016). Inconsistencies of the above parameters calculated for the corresponding FHR\_U and FHR\_E signals were estimated using the symmetric mean percentage difference SMPD, where the differences between values are related to their mean. Since the normality assumption was verified using the Shapiro-Wilk test, the statistical significance (using paired Student's t-test) of the differences between the corresponding parameters obtained in the FHR\_U and FHR\_E signals was examined.

### Indirect Signals Comparison via Beat-to-Beat Variability

It is generally believed that a very high predictive value in relation to the early detection of the fetal distress is provided by the instantaneous FHR variability parameters (Cesarelli et al., 2009). They are determined from the FHR signal in a form of the time event series—a sequence of events unevenly located in time, providing the successive cardiac cycles duration expressed in milliseconds.

This format significantly differs from that available at the output of a fetal monitor—the FHR values evenly spaced at every 250 ms. This measurement period has been established to be not longer than the shortest physiologically allowed heart cycle, however with characteristic information redundancy for low FHR values (e.g., the FHR value equal to 50 bpm is represented by four duplicated subsequent values; Lee et al., 2009; Goncalves et al., 2013).

Therefore, the FHR\_U and FHR\_E signals were subjected to reconstruction of the above mentioned time event series representation. This procedure relied on taking from the evenly distributed time series, the values according to the timing signal being constituted by the fetal QRS complexes additionally obtained from the KOMPOREL System. The resulting signal is a sequence of time-ordered events corresponding to subsequent occurrences of the fetal QRS complexes (or more precisely the R-waves). Created according to the described procedure the FHR\_U and FHR\_E signals in the form of time event series, provided the basis values for determination of the instantaneous variability indices. These indices, widely acclaimed in the literature (Romano et al., 2016b), quantitatively describe the long- and short-term FHR variability.

In this study the following indices were analyzed: Haan\_LTI, Haan\_STI, Yeh\_II, Yeh\_STI, Organ\_LTV, Organ\_STV, Dalton\_LTV, Dalton\_STV, Zugaib\_LTV, Zugaib\_STV. In order to standardize the results the indices were determined within 1-min segments (Kubo et al., 1987; Jezewski et al., 2006). Each consecutive segment comprised only those instantaneous FHR values which were determined using the heart beats contained in the given segment (Cesarelli et al., 2009).

While processing a given segment, if percentage of valid values, relevant for a given index (according to its definition), was less than 20%, it was assumed that the index value was undetermined for that minute. Such cases are the result of the signal loss episodes in the analyzed signals. In a series of minute values calculated for a given signal, they are defined as a 1-min loss of the given index value and marked with the value of −1. The index average value for a given signal is calculated from all 1-min values, excluding those marked as undetermined. Additionally, a non-linear parameter of instantaneous FHR variability was proposed, in a form of the regularity measure—the sample entropy index (SampEn) (Signorini and Magenes, 2014). It was determined in the windows covering 300 heart events (FHR values), and expressed in milliseconds as a measure of period. The parameters of SampEn function were set at: dim = 1 and r = 0.1. For the given signal, the SampEn index represents the mean value of sample entropy determined in successive windows. Inconsistencies of indices describing the instantaneous FHR variability, determined for corresponding FHR\_U and FHR\_E signals were evaluated with the SMPD, where the difference between values are related to their mean. The statistical significance (using paired Student's t-test) of the differences between the values of the indices obtained for each signal pair was also examined.

### Indirect Signals Comparison via Fetal Outcome Prediction

Analysis of the FHR signal leads to its classification as corresponding to normal or abnormal fetal state. Since at the time of fetal monitoring, there is no other diagnostic method which could be able to confirm a correctness of the signal classification, the FHR signals, being acquired during pregnancy, are retrospectively assigned to true fetal outcome (newborn state) (Chudacek et al., 2014; Romano et al., 2016a). It is justified, as in obstetrics it is assumed that the normal fetal outcome has to be result of proper fetal development during the pregnancy period. Excluding the cases when the labor process itself caused negative effects for the fetal state, the same assumption can be applied for abnormal fetal outcome. It is generally believed that this relationship is maintained in case of deliveries by cesarean section due to the maternal reasons. In the collected database the vast majority of cases were of this type. The 60 patients (recordings) were classified as belonging to a normal or abnormal group using the information on fetal outcome. The abnormal state was set if at least one of the following conditions was met: Apgar score (at 5 min) <7, pH <7.2, BE >12, NICU stay > 24 h, or birth weight percentile <5%. Finally, the research material included 19 recordings with abnormal state and 41 with normal state assigned retrospectively. It is important that for the majority of patients with abnormal fetal outcome, the pregnancy was ended by cesarean section due to the maternal indications (16 cases). It suggests that course of delivery imposed no negative effect on the newborn. The database contains six recordings for which the abnormal state was set due to three or more conditions met.

The analysis of ability for prediction of the fetal outcome was performed separately for the ultrasound and the fetal electrocardiography approaches. Each FHR\_U (FHR\_E) signal was represented by the feature set comprising: 15 parameters determined by the MONAKO System, 10 indices describing the instantaneous FHR variability and SampEn entropy measure. The differences between the values of these features obtained in two groups (normal and abnormal fetal outcome) were expressed by the mean percentage difference (MPD) for both groups where the normal outcome group was taken as reference. The statistical significance of the difference between the mean values of each feature obtained for two groups was assessed using Student's t-test.

The indirect comparison of FHR\_U and FHR\_E signals, as for predicting the fetal outcome, was based on the capability to classify the FHR signals from a given acquisition methods into normal and abnormal analyzing the determined clinical FHR parameters.

#### RESULTS

Comparison analysis considering two different methods of FHR signal acquisition has been carried out using 60 pairs of FHR\_E and FHR\_U signals, which were obtained during the monitoring sessions of 60 patients. After a full signals synchronization and trimming, the total length of recordings was equal to 1995 min, with an average length of 33.3 min (SD = 10.9 min). The recordings were characterized by low signal loss (details in **Table 1**). For the US method the mean signal loss was equal to 4.5%, whereas for the FHR signals obtained via FECG, the loss level was more than two times lower, which is expressed by the mean signal loss of only 1.8%. High signal loss was observed in few recordings, especially for FHR\_U. It could be noted that for more than half of the recordings, the signal loss for both methods did not exceed 1%. Such low level of the signal loss and the sufficient length of individual recordings enabled further estimation of the inconsistency between both acquisition methods, using the mean values of clinical quantitative parameters, determined during the FHR signal analysis.

A direct comparison of the FHR\_E and FHR\_U signals has been based on estimation of the differences between the corresponding instantaneous FHR values (provided every 250 ms). It enabled the metrological assessment of the inconsistency between the two acquisition methods. Comparison at that stage was performed for each individual recording and the summary of descriptive statistics for entire research material is presented in **Table 1**. The mean difference value (MD) obtained for all recordings was −0.23 bpm, which in relation to an average value of FHR (about 140 bpm) gives a relative error of 0.2%. It means that the measurement bias between the two methods does not occur. When analyzing the MD values calculated for particular recordings we noted that it did not depend on the measured FHR signal value. It is shown by the Bland-Altman plot of MD values against the average FHR values (from FHR measured at 250 ms) obtained for particular recordings (**Figure 2**). From the point of view of both visual and automated assessment of the FHR signal variability, the more important seems to be the MAD, which was equal to 1.24 bpm for all the considered signals. Such value does not affect a visual evaluation of the signal, since it is lower than the printing resolution of the FHR waveforms, as well as the resolution of a human eye.

In the computer-aided system, the automated analysis aimed at determination of clinically important FHR patterns is carried out using the FHR values averaged over 2.5 s. It makes the comparison between FHR\_E and FHR\_U signals represented by such averaged values especially important. Averaging processes caused a slight decrease of the signal loss in both types of

TABLE 1 | Descriptive statistics of the signal loss and values of differences MD, MAD, determined for the final set of recordings from the electrocardiography (FECG) and the ultrasound method (US), where the FHR signals were expressed as the original 250 ms measures, and as the values averaged over 2.5 s periods.


signals, to the values: 4.3% in US, and 1.4% in FECG (**Table 1**). Comparison at this stage was performed for particular recordings and the summary of descriptive statistics for entire research material is presented in **Table 1**. In general, inconsistency between FHR signals after averaging decreased. In this case, the bias between the two methods also does not occur. The MAD value has decreased significantly, to only 0.7 bpm, which is below the 1 bpm level—the minimum accuracy of the FHR measurement. In relation to the average value of the FHR signal (140 bpm), the relative inconsistency was equal to 0.5%. Considering these results we could assume that such small inconsistency should lead to the similar values of clinical parameters provided by both methods. However, we have to keep in mind that some of these parameters are particularly sensitive to the FHR changes, being a result not only of the mean FHR difference, but rather of the distribution of FHR differences in time. In contrast, some other parameters of the FHR signal are sensitive to temporary high differences in the FHR signals. So, in case of these parameters the differences between both methods may occur. Such formulated assumptions have been verified in the next stage of the inconsistency analysis the indirect comparison of both signals. It was based on the interpretation of the differences between the particular FHR signal parameters, which are provided by an automated analysis in the fetal monitoring system.

Descriptive statistics (mean values and SD) for individual parameters describing quantitatively the FHR signal, which have been determined for signals from both methods, are presented in **Table 2**. For each parameter the differences between the FHR\_E and FHR\_U signals, were assessed using the SMPD. It was justified because for any of those parameters no significant difference between them was noted. As it could be expected, the low SMPD value of 0.1% obtained for M\_FHR and M\_BL parameters (being significantly dependent on an averaging process of the FHR measurements), was similar to the relative inconsistency reported in the direct signals comparison. In turn, the SMPD values calculated for the following parameters: HE\_D, LE\_D, STV, ACC, DEC, OSC\_IV, differ significantly from the values reported for other parameters, and even more TABLE 2 | Values of clinically important parameters of quantitative description of FHR-E and FHR-U signals, from FECG and US, obtained in a computer-aided fetal monitoring system, together with the symmetric mean percentage difference SMPD estimating the inconsistencies between both the methods.


\**p* < *0.05 (paired t-test).*

from the results of the direct comparison. It confirms the above mentioned assumption that some parameters (e.g., STV) are sensitive to distribution of the FHR differences in time, whereas another ones (e.g., DEC) to a temporary high difference value. Particularly, high SMPD = 54% noted for the number of recognized decelerations, is a results of two factors: the direct differences between the signal values and the differences in the signal loss episodes. The signal loss for FHR\_U is on average twice higher than for FHR\_E. In addition, the autocorrelation technique, commonly used in the US method to determine the signal periodicity, is often not able to follow the rapid decrease of FHR signal related to deceleration, which results in signal loss episodes (**Figure 3**). This, in turn, causes that the deceleration is not recognized, because it does not meet the established criteria of amplitude and duration.

Indirect comparison of the FHR\_U and FHR\_E signals was performed on the basis of the variability indices defined for signal represented as time event series—the heart beats. Summary of the results (mean values, SD, and SMPD) for the selected 11 parameters describing the FHR signal variability is presented in **Table 3**. The results clearly show that the FHR\_E signal is characterized by higher variability then the FHR\_U. Inconsistencies for the long-term variability indices were at about 10%, whereas for the short-term indices they were five times higher, reaching about 50%. For the short-term variability we noticed significant difference for both methods.

TABLE 3 | Results of the FHR\_E and FHR\_U signal analysis, concerning the long- and short-term variability, calculated using signal in a form of time event series—a sequence of events unevenly localized in time, together with the SMPD values estimating the inconsistencies between both the methods.


\**p* < *0.03;* <sup>+</sup>*p* < *0.001;* #*p* < *0.0001 (paired t-test).*

With regard to such large inconsistencies it has to be decided which of the two methods may be considered as providing the FHR variability description being closer to the true one. The answer is not obvious, because in this work no reference signal was acquired simultaneously with two analyzed methods. Such gold standard can be provided by previously mentioned the direct fetal electrocardiography, where the pure FECG is acquired from the fetal head. In the previous studies where the FHR signal from ultrasound method was compared with the reference one, it has been shown that ultrasound method underestimates the short-term variability on a level between 20 and 40% in reference to direct fetal electrocardiography. A similar trend can be seen in **Table 3**, where the FHR\_U is compared with FHR\_E obtained from the abdominal fetal electrocardiogram.

The signal database has been divided into two groups: normal and abnormal, according to the established fetal outcome criteria. The description of two signal groups is shown in **Table 4**. Summary of the FHR signal analysis results, comprising 15 clinical parameters determined for both groups of fetal outcome, are shown separately for the ultrasound method (**Table 5**) and abdominal electrocardiography (**Table 5**). In addition to the mean value and standard deviation, the mean percentage difference MPD was calculated, assuming the values obtained in normal fetal outcome group as the reference. Apart from assessing the statistical significance of the difference between the normal and abnormal groups, also the tendency of changes was studied for the particular parameters between these groups. It was carried out to check whether the observed tendency would be consistent with the clinical interpretation of those parameters. Considering the number of ACC patterns, a significant difference between the abnormal and normal groups was noted, and the tendency was consistent. But for the number of DEC patterns no significant difference was observed. The MPD took the opposite values: negative value for ultrasound (−48%) and positive for electrocardiography (38%). Thus, a clinical interpretation of decelerations as a sign of fetal distress has been confirmed

FIGURE 3 | An example of a 12-min fragment of signal pair—result of an indirect comparison to evaluate the impact of the FHR signal acquisition method on clinically relevant parameters, determined by a computer-aided fetal monitoring system. The autocorrelation technique, commonly used in the US method to determine the signal periodicity, is often not able to follow the rapid decrease of FHR\_U signal related to deceleration, which results in signal loss episodes. This, in turn, causes that the deceleration is not recognized, because it does not meet the established criteria of amplitude and duration. Graphic markers of the analysis results illustrate the signal loss (above the curve), the estimated FHR baseline (line fitted on FHR curve) and detected deceleration episodes (horizontal bars under the curve).

TABLE 4 | Descriptive statistics of the selected parameters, describing the labor and fetal outcome, presented separately for two groups of the research material, including 19 recordings with abnormal and 41 recordings with normal fetal outcome.


only in case of electrocardiography. It is mainly caused by the different amount of the signal loss episodes in both types of signal. Particularly, in the FHR\_U the signal loss often occurs during the deceleration which leads to misdetection of this pattern (**Figure 3**). As for the accelerations both methods follow the clinical meaning as they provided a higher number of such episodes relating to fetal well-being in normal group than in abnormal one.

As above, no statistically significant differences were noted for all the indices describing both the short-term and long-term FHR variability. Decrease of most of the long-term variability indices was noted in abnormal groups (MPD ranged from −20 to −9%) for both FHR\_U (**Table 6**) and FHR\_E signals (**Table 6**). This is consistent with the clinical interpretation of these indices, since a decrease of FHR variability is regarded as the fetal distress sign. Also, for most of the short-term variability indices a decrease was noted in the abnormal group. In the case of the US method, the MPD takes value of −12%, but for the FECG method it ranges between smaller values—from −4 to −2%. An interesting property is shown by the Geijn\_STV index, which regardless of the acquisition method used, exceeds the mentioned range for short-term variability indices in the abnormal group, as it is significantly reduced (MPD = −36% for both methods). Contrary to the Geijn\_STV, the Haan\_STI although also exceeds the mentioned range for FHR\_E signals, but shows different tendency—it is overestimated by 3%. Finally, it should be noted that the variability indices are capable to differentiate the signals relating to normal and abnormal fetal state, and among them the Geijn\_STV index is particularly effective.

#### CONCLUSIONS

The paper proposes an extended process of comparing two different methods of fetal heart rate monitoring. It takes into account the issues associated with comparing different biomedical signals—the unsatisfying usability of the results obtained from direct signal comparison. High reliability of the comparison results can only be ensured by using the clinically significant parameters determined for signals acquired by both methods being analyzed (Jezewski et al., 2006). Initial preparation of the research material has been stated as very important step too, in order to ensure the results are not affected with different measurement conditions. Two different methods of measuring the fetal heart rate signal were analyzed: indirect electrocardiography for recording the electrical heart activity from maternal abdominal wall, and the pulsed Doppler ultrasound method based on mechanical activity of the fetal heart (Cohen et al., 2012). None of these methods can be considered as a reference, due to a number of measurement error sources identified (Goncalves et al., 2015).

There is a strong conviction that the ultrasound method, as leading in clinical practice, can serve as a quasi-reference. It applies when the acquired fetal heart rate signal is interpreted both visually and more and more often, using the results of quantitative analysis in the computer-aided fetal monitoring system. On the other hand, we observe growing expectations on the indirect fetal electrocardiography, especially with regard to its application for pregnancy telemonitoring at home (Martinek et al., 2015). The above issues justified a need for the comparison between electrocardiography and ultrasound method presented in this paper in a view of usability of the results obtained. The complex process of data synchronization and validation within the research material resulted in 60 pairs of FHR signals, with an average duration of about 33 min. Obtained low level of the signal loss (4.5% for the US and 1.8% for FECG method) enabled to perform both direct comparison and indirect one—by using clinically relevant parameters (Reinhard et al., 2010; Wrobel et al., 2015b).

From direct comparison it has been resulted there is no measurement bias between the acquisition methods. The mean difference between the FHR\_E and FHR\_U signal was equal to −0.2 bpm (SD = 0.38 bpm). On the other hand, the mean absolute difference measured between the methods, being important for both visual and computer-aided signal analysis, was equal to 1.2 bpm. When relating to typical FHR level of about 140 bpm, the inconsistency takes the value of about 0.9%. These results are similar to those obtained in Cohen et al. (2012) where the ultrasound method was compared with the reference direct fetal electrocardiography and in Jezewski et al. (2012), where abdominal electrocardiography was compared with reference in a similar manner.

Such low differences do not affect the visual assessment of the FHR signal, taking into account a resolution limit of both a printer and human eye (Reinhard et al., 2012). However, when analyzing the results of the indirect comparison, by using the parameters quantitatively describing the clinical features, the inconsistencies of several percent were noted. This particularly affects the patterns being sensitive to the instantaneous differences in FHR values, like acceleration (7.8%) and particularly deceleration (54%), where (for the ultrasound method) the signal loss within the episodes has a significant impact (Voicu et al., 2010). Similarly, significant differences were


TABLE 5 | Summary of the FHR-E and FHR-U signal analysis comprising 15 clinical parameters determined using computer-aided fetal monitoring system, for two groups of fetal outcome, together with the mean percentage difference MPD depicting the inconsistencies between both the groups.

\**p* < *0.05;* #*p* < *0.05 (t-test).*

TABLE 6 | Results of the FHR\_E and FHR\_U signal analysis, in terms of the instantaneous FHR variability assessment, calculated for signal in a form of time event series, together with the mean percentage difference MPD depicting the inconsistencies between both the groups.


noted between the ultrasound method and the reference direct electrocardiography in Desai et al. (2013).

The results obtained for the long- and short-term FHR variability indices show significant overestimation of their values in the signal acquired using the electrocardiography, by 10% and 50%, respectively. However, the results of work (Jezewski et al., 2006), where the US method was related to the reference direct FECG, have shown that the electrical method provides significantly higher FHR variability. Hence it leads to conclusion, that the variability indices for the FECG signal acquired from the abdominal wall represent the true FHR variability, whereas in the ultrasound method they are significantly underestimated (Goncalves et al., 2013). On the other hand, a comparison of these methods, through a clinical interpretation of the FHR signals for fetal outcome prediction was examined. It showed that ability of various clinical parameters to distinguish between normal and abnormal groups do not depend on the acquisition method. That was confirmed by similar tendency of changes of clinical parameters determined in both groups.

In summary we can conclude that the abdominal fetal electrocardiography, being considered as an alternative to the ultrasound based approach for certain application, like a pregnancy telemonitoring at home, does not change significantly the interpretation of the FHR signal. Equivalence of these methods was confirmed for both visual assessment and automated analysis of the signals. Despite the lack of reference signal, it can be proved indirectly that the abdominal fetal electrocardiography provides more reliable description of the instantaneous FHR variability. It's another advantage over the ultrasound method relates to a lower signal loss. However, this conclusion coming from analysis of the signals collected in hospital conditions may undergo final verification when both methods will become widely applied in the systems for pregnancy telemonitoring.

#### ETHICS STATEMENT

The study protocol was approved by the Ethical Committee of the Silesian Medical University, Katowice, Poland (NN-013-123/03). Subjects read approved consent form

#### REFERENCES


and gave written informed consent to participate in the study.

### AUTHOR CONTRIBUTIONS

JJ, JW, and AM designed the study approach and experiments. KH and AM described the research background in the introduction section. KH and JJ collected the research material. TK and RM provided the signal synchronization procedures. JJ, JW, and MJ were responsible for analysis and interpretation of data. TK and JJ performed statistics. JW wrote the manuscript with contributions from JJ, KH, RM, and MJ. All authors read and approved the final manuscript.

#### ACKNOWLEDGMENTS

This scientific research work is supported by The National Centre for Research and Development of Poland (grant No STRATEGMED2/269343/18/NCBR/2016).


recordings. IEEE T Bio. Med. Eng. 63, 1269–279. doi: 10.1109/TBME.2015. 2493726


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer HA and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.

Copyright © 2017 Jezewski, Wrobel, Matonia, Horoba, Martinek, Kupka and Jezewski. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Assessment of Fetal Development Using Cardiac Valve Intervals

Faezeh Marzbanrad<sup>1</sup> \*, Ahsan H. Khandoker 2, 3, Yoshitaka Kimura<sup>4</sup> , Marimuthu Palaniswami <sup>2</sup> and Gari D. Clifford5, 6

<sup>1</sup> Department of Electrical and Computer Systems Engineering, Monash University, Clayton, VIC, Australia, <sup>2</sup> Electrical and Electronic Engineering Department, University of Melbourne, Melbourne, VIC, Australia, <sup>3</sup> Biomedical Engineering Department, Khalifa University of Science, Technology and Research, Abu Dhabi, United Arab Emirates, <sup>4</sup> Graduate School of Medicine, Tohoku University, Sendai, Japan, <sup>5</sup> Department of Biomedical Informatics, Emory University, Atlanta, GA, United States, <sup>6</sup> Department of Biomedical Engineering, Georgia Institute of Technology, Atlanta, GA, United States

Edited by:

Zbigniew R. Struzik, University of Tokyo, Japan

#### Reviewed by:

João Francisco Bernardes, University of Porto, Portugal Linda E. May, East Carolina University, United States

> \*Correspondence: Faezeh Marzbanrad faezeh.marzbanrad@monash.edu

#### Specialty section:

This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology

Received: 29 January 2017 Accepted: 01 May 2017 Published: 17 May 2017

#### Citation:

Marzbanrad F, Khandoker AH, Kimura Y, Palaniswami M and Clifford GD (2017) Assessment of Fetal Development Using Cardiac Valve Intervals. Front. Physiol. 8:313. doi: 10.3389/fphys.2017.00313 An automated method to assess the fetal physiological development is introduced which uses the component intervals between fetal cardiac valve timings and the Q-wave of fetal electrocardiogram (fECG). These intervals were estimated automatically from one-dimensional Doppler Ultrasound and noninvasive fECG. We hypothesize that the fetal growth can be estimated by the cardiac valve intervals. This hypothesis was evaluated by modeling the fetal development using the cardiac intervals and validating against the gold standard gestational age identified by Crown-Rump Length (CRL). Among the intervals, electromechanical delay time, isovolumic contraction time, ventricular filling time and their interactions were selected in a stepwise regression process that used gestational age as the target in a cohort of 57 fetuses. Compared with the gold standard age, the newly proposed regression model resulted in a mean absolute error of 3.8 weeks for all recordings and 2.7 weeks after excluding the low quality recordings. Since Fetal Heart Rate Variability (FHRV) has been proposed in the literature for assessing the fetal development, we compared the performance of gestational age estimation by our new valve-interval based method, vs. FHRV, while assuming the CRL as the gold standard. The valve interval-based method outperformed both the model based on FHRV. Results of evaluation for 30 abnormal cases showed that the new method is less affected by arrhythmias such as tachycardia and bradycardia compared to FHRV, however certain types of heart anomalies cause large errors (more than 10 weeks) with respect to the CRL-based gold standard age. Therefore, discrepancies between the regression based estimation and CRL age estimation could indicate the abnormalities. The cardiac valve intervals have been known to reflect the autonomic function. Therefore the new method potentially provides a novel approach for assessing the development of fetal autonomic nervous system, which may be growth curve independent.

Keywords: fetal development, gestational age, 1D Doppler ultrasound, cardiotocography (CTG), fetal electrocardiography (fECG), autonomic nervous system (ANS), fetal monitoring, systolic and diastolic time intervals

# 1. INTRODUCTION

Estimation of the Gestational Age (GA) is crucial for antenatal diagnosis, monitoring fetal growth and detecting Intra-Uterine Growth Retardation (IUGR), predicting the delivery date and management of pre-term and post-term pregnancies, to ultimately prevent perinatal and neonatal mortality (Alexander et al., 1996; Taipale and Hiilesmaa, 2001; Bhutta et al., 2014; Chauhan et al., 2014). It is also a fundamental factor in ensuring the safety and effectiveness of medications during pregnancy (Reis and Källén, 2010; Andersen et al., 2013; Li et al., 2013). The GA has been traditionally estimated based on the Last Menstrual Period (LMP), which is the most affordable method specially for low and middle income countries (Wang et al., 2011; Deputy et al., 2017). However it is subject to human errors in recall or data entry, as well as biologically associated errors (Dietz et al., 2007; Lynch and Zhang, 2007). For example, the assumption of a regular 28-day menstrual cycle and ovulation on 14 days after the first day of LMP are not consistent and may vary from case to case (Dietz et al., 2007). As reported in the literature, the clinically estimated GA as collected on certificates of live birth based on prenatal and neonatal clinical assessments, exceeds the LMP-based GA by 2 weeks or more for more than 40% of the cases (Alexander et al., 1995). A recent study found the ovulation day as the most accurate predictor compared to LMP and ultrasound-based methods (Mahendru et al., 2016).

A more accurate and reliable estimation of the GA is provided through obstetric ultrasound which has been clinically established as the gold standard (Lynch and Zhang, 2007; Papageorghiou et al., 2014). A variety of sonographic measurements including Biparietal Diameter (BPD), Crown-Rump Length (CRL), Head Circumference (HC), Abdominal Circumference (AC) and Femur Length (FL) are used to estimate the GA (Hadlock et al., 1982; Dietz et al., 2007; Lynch and Zhang, 2007; Papageorghiou et al., 2014). Although these measures provide a more reliable estimation of GA compared to LMP, they are all based on physical growth (mass or proportions), which is affected by genetic variations (e.g., head size and shape in fetuses), gender and inherent variability in the fetal growth process (Hadlock et al., 1981; Sherwood et al., 2000; Lynch and Zhang, 2007; Kullinger et al., 2016). These methods may also systematically overestimate or underestimate the GA of the fetuses which are respectively large or small for GA (Sherwood et al., 2000; Lynch and Zhang, 2007). Unsuitable positioning of the fetus during measurement also causes error and the technique is subject to operator error, and the quality of the images (Hunter, 2009; Callen, 2011). For example, 95% confidence intervals of ±4 weeks were found for FL, which is one of the most accurate estimators. The prediction interval may be as large as ±7 weeks for the estimators with higher standard errors, such as AC (Sherwood et al., 2000). The error also increases with the gestational age and generally the ultrasound methods are more precise when performed in the first-trimester (Caughey et al., 2008; Falatah et al., 2014; Al-Amin et al., 2015). Pathological conditions may also introduce a high levels of inaccuracy or significant bias in many estimation methods.

Although in high income countries routine skilled ultrasound screening is performed, many factors limit its use in low income countries, including high cost of the equipment, lack of trained sonographers or physicians, as well as the skill required to perform a GA estimation test (Wang et al., 2011; McClure et al., 2014). We therefore propose an alternative technique to be used as an adjunct in estimating the GA where ultrasound imaging methods are unavailable or inadequate due to pathologies, unsuitable positioning, limited skills and technical issues.

One promising alternative GA estimator is Fetal Heart Rate (FHR) (Cha et al., 2001; Hoyer et al., 2013; Tetschke et al., 2016). FHR can be measured with affordable apparatus and little need for prior skill, and hence is a feasible approach in low income countries (Tezuka et al., 1998; Stroux et al., 2014). Early studies found a comparable accuracy of FHR-based method with CRL, in early pregnancy (38–64 days) (Tezuka et al., 1998), and the methods were further improved recently being more focused on neurological development (Hoyer et al., 2013; Tetschke et al., 2016). While ultrasound-based techniques are generally based on the physical development and influenced by genetic variations, FHR provides a marker for neuro-physiological development of the fetus, since it reflects the Autonomic Nervous System (ANS) control of the cardiovascular system, which matures through the progress in pregnancy. Various linear, nonlinear and complexitybased FHR variability (FHRV) parameters have been found to be closely related to the fetal development (Van Leeuwen et al., 2003; Hoyer et al., 2009; Wallwitz et al., 2012; Hoyer et al., 2013; Tetschke et al., 2016). Vagal and sympathetic activity rhythms and their interactions has been traditionally attributed to different frequencies of FHR fluctuations and their ratios (European Society of Cardiology, 1996) which can be evaluated during fetal development. More recently, a "functional Fetal Autonomic Brain Age Score" (fABAS) was introduced which leverages the FHR patterns in a multivariate analysis (Hoyer et al., 2013). However, FHR is influenced by arrhythmias, fetal behavioral/sleep states and heart rate patterns such as FHR accelerations and even maternal psychological and physiological conditions, particularly in mid- and late-gestation (Mantel et al., 1991; Monk et al., 2000; Ivanov et al., 2009; Marzbanrad et al., 2015b). These factors complicate the assessment of fetal development based on FHR.

Fetal cardiac valve intervals are alternative measures which could be obtained from non-invasive, low cost and easy-to-operate devices, and used as reliable markers for fetal development and well-being (Shakespeare et al., 2001; Khandoker et al., 2009; Marzbanrad et al., 2013b). These intervals are based on the opening and closing time of the fetal cardiac valves, namely the atrioventricular and semilunar valves. Automated techniques for estimation of these intervals from non-invasively recorded one-dimensional Doppler Ultrasound (1-D DUS) signal (conventionally used as FHR monitor) were proposed in our previous papers (Marzbanrad et al., 2013b, 2014a). The valve intervals can also be used to assess the ANS function, as an alternative to the FHRV, since the cardiac mechanics are known to reflect the autonomic control in the literature on adults (Berntson et al., 1994; Cacioppo et al., 1994a,b; Di Rienzo et al., 2013). As an example, the Pre-ejection Period (PEP), which is the interval from the onset of ventricular depolarization (the beginning of the QRS complex on the electrocardiogram) to the opening of aorta, reflects sympathetic influences on the heart (Cacioppo et al., 1994a; Mensah-Brown et al., 2010). In previous studies we found significant changes in the valve intervals with advancing GA (Marzbanrad et al., 2013a,b). In the work presented here we hypothesize that the fetal cardiac valve intervals, which are estimated automatically, can be used as a novel alternative measure of the GA, reflecting the physiological development of fetus.

# 2. METHODS

#### 2.1. Subjects and Data Acquisition

Doppler ultrasound and abdominal ECG signals were recorded simultaneously at Tohoku University Hospital, Japan, from 57 pregnant women with healthy single pregnancy who were not under any medication and 30 cases with fetal arrhythmia or abnormalities. The type of abnormalities and the GA for these 30 cases are presented in **Table 4** and more details about these arrhythmia and abnormalities can be found in Murray (2007), Allan et al. (2000), Allen et al. (2013), and Abuhamad and Chaoui (2012). All 87 fetuses had a GA of between 16 and 41 (32 ± 6) weeks at the time of recording. The GA was estimated using ultrasound imaging by a trained sonographer, by measuring fetal CRL at about 10 weeks, which is the length of embryos and fetuses from the head top (crown) to the bottom of the buttocks (rump). The study protocol was approved by Tohoku University Institutional Review Board and written informed consent was obtained from all participants. **Table 1** summarizes the CRL-based gestational age, maternal age, weight and height for the healthy and abnormal cases.

The 1-D DUS signal was generated using a 1.5 MHz Corometrics 5700 Ultrasound transducer and the abdominal ECG signals were collected by a multichannel data acquisition system (fetal monitor 116, Corometrics Medical Systems Inc) with 1,000 Hz sampling frequency and 16 bit resolution. Twelve electrodes were used for abdominal ECG recordings, ten of which were arranged on the mother's abdomen, one reference electrode on the back and one electrode was set at the right thoracic position. The DUS transducer was placed on the lower abdomen and the audio output was connected to the input channel of the fetal monitor. All DUS and ECG recordings were 1 min in length and sampled at 1 kHz with 16-bit resolution. More details about the experimental set up can be found in Sato et al. (2007).

TABLE 1 | Maternal age (years), height (cm) and weight (kg) as well as the CRL-based GA (weeks) for normal and abnormal groups are presented as mean ± standard deviation.


## 2.2. fECG Extraction

Data from 12 channels were recorded bipolarly from the electrodes placed on the maternal abdomen, sampled every 1 ms (1 kHz sampling) with 16-bit resolution and bandpass filtered by 1–100 Hz finite impulse response filter. To separate fECG from the composite abdominal signal, a combination of maternal ECG cancellation and Blind Source Separation with a Reference (BSSR) was employed (Sato et al., 2007). In brief, electrical activities of the heart can be modeled as a vector in the direction of excitation, which is sometimes called the heart vector (Symonds et al., 2001). The maternal ECG component was excluded by subtracting the linear combination of mutually orthogonal projections of the heart vector. Subsequently, BSSR was used to extract fECG from complex mixture using DUS signal as a reference (Sato et al., 2007). Fetal QRS locations were detected by a modified Pan and Tompkins peak detection method as described in Behar et al. (2013).

# 2.3. Evaluation of Data Quality

The quality of 1-D DUS and fECG signals were assessed to exclude low quality signals and to evaluate the relationship of the final error in the GA estimation with the quality scores.

#### 2.3.1. fECG Signal Quality

A state of the art Signal Quality Index (SQI), known as "bSQI," was used. This metric evaluates the agreement between two QRS detection methods with different robustness to noise (Clifford et al., 2012). The bSQI metric takes a range between 0 (lowest quality) and 1 (highest quality).

#### 2.3.2. 1-D DUS Signal Quality

Quality assessment of the DUS signal was performed using a method described in our earlier work (Marzbanrad et al., 2015a). The method is based on various quality indices of the high frequency component of the DUS signal. To isolate the high frequency component which is linked to the valves' movements, the DUS signal was decomposed by continuous wavelet analysis, as described in our earlier work (Khandoker et al., 2009; Marzbanrad et al., 2014a). Using a second order complex Gaussian as the mother wavelet, the signal at scale 2 (∼ 200 Hz) was extracted and smoothed. The envelope of the absolute value of this signal was then estimated by interpolating the maxima and smoothing with a low pass filter. Each envelope was segmented into cardiac cycles using the corresponding RR intervals, estimated from fECG. The signal segments were then normalized by subtracting the mean and dividing by the standard deviation.

Twelve features were selected mainly based on the signal properties in the valve motion ranges compared to the remaining time intervals. The plausible valve motion ranges were defined as: Mc: (9–44), Ao: (45–90), Ac: (200–260), Mo: (265–326), all in msec following the segment onset (the preceding R-peak) (Khandoker et al., 2009; Marzbanrad et al., 2013b). The features selected were as follows and all were normalized to [0 − 1]:

• The ratio of the power (SQI1), number of peaks (SQI2), mean peak amplitude (SQI3) and variance (SQI4) in the valve motion range to the values in the remaining time intervals.


An overall quality metric was obtained from the features SQI1,2,..,12 using a Naïve Bayes (NB) classifier with kernel density estimate. The classifier was trained based on 345 cardiac cycles of the DUS signals which were annotated for quality by four independent annotators, as described in our previous paper (Marzbanrad et al., 2015a). The NB classifier uses the training data to estimate the conditional distribution of the features given the classes and also distribution of the classes. Then it estimates the posterior probability through the Bayes rule and classifies each sample to the most probable class. The same trained classifier was used in this study to classify the DUS quality. In our previous paper we used 10-fold cross validation and found the accuracy of 0.86 and 0.84 in train and test set, respectively.

## 2.4. Estimation of Cardiac Valve Intervals

The cardiac valve intervals are illustrated in **Figure 1**. These intervals were obtained based on the onset of the ORS complex detected as described above, and the opening and closing of the valves detected from the high frequency component of the DUS signal. The valve motion events were detected using a model-based method that was presented in our previous work (Marzbanrad et al., 2014a). This latter method is now summarized. The envelope of the high frequency component of the DUS signal, which were normalized and segmented into cardiac cycles (as described in Section 2.3.2), were clustered into six different patterns using K-means clustering. The key idea was to find the following events which correspond to the peaks of the high frequency components: Aortic valve opening (Ao), transitional event (T1), Aortic valve closing (Ac), transitional event (T2), Mitral opening (Mo), transitional event (T3), Mitral closing (Mc), transitional event (T4). The transitional events are related to extra peaks that do not correspond to any valve motion. A hybrid Support Vector Machine-Hidden Markov Model (SVM-HMM) was trained for each cluster separately, using the time (phase) and amplitude of the peaks of the signal as features, corresponding to one of the valve motion or transition events. The training and validation of this approach were based on expert annotation and simultaneous fetal echocardiography images, and were carried out in our earlier work (Marzbanrad et al., 2014a). To identify the events, each segment of the normalized envelope of the high frequency component was matched to the clusters that for which it had the minimum Euclidean distance to cluster's centroid. Then the sequence of events, which were attributed to the peaks of the signal, were identified by the Viterbi algorithm using the trained SVM-HMM specific to the corresponding cluster. The block diagram of this method is shown in **Figure 2** and more details can be found in Marzbanrad et al. (2014a).

# 2.5. Estimation of the Gestational Age

Three sets of parameters were used to estimate the GA:

• Valve-timing parameters: From the parameters shown in **Figure 1**, Electromechanical Delay Time (EDT), Isovolumic Contraction Time(ICT), Ventricular Ejection Time (VET), Isovolumic Relaxation Time (IRT) and Ventricular Filling Time (VFT) were selected. Only Pre-Ejection Period (PEP)

and Systolic Time Interval (STI) were excluded as they were linearly related to other intervals.


In order to estimate the GA from these parameters a stepwise regression analysis was employed based on individual and all combinations of parameters and the models including an intercept, linear, squared terms and cross-products. Stepwise regression automatically adds to or removes from the model in a forward and backward process to determine a final model, using an F-test applied to the sum of the squared error before and after adding a parameter (p < 0.05) as the criterion for including a parameter. Root Mean Squared Error, R-squared, adjusted Rsquared and the F-test results vs. constant model were calculated for the regression of each set of parameters. An average leaveone-out cross-validation error in GA estimation was calculated. The difference between the CRL-based and regression-based GA estimate was made at every stage to provide an estimate of out of sample performance of the proposed approach.

The GA estimation error was compared for different parameters, including: FHR, fECG quality score and DUS signal quality score. The improvement of the GA estimation by applying threshold on the quality of the signals was also evaluated. Finally the optimal regression model which was obtained based on healthy cases was then used to estimate the GA of abnormal cases.

#### 3. RESULTS

#### 3.1. Stepwise regression results

Using fetal heart valve timings as parameters, stepwise regression resulted in the following regression model for all healthy fetuses, without excluding the cases with low quality signals:

$$\text{Estimated GA} = a\_0 + a\_1EDTA + a\_2ICT + a\_3VFT + a\_4EDTA$$

$$\astICT + a\_5ICT \ast VFT$$

where a0, a1, ..., a<sup>5</sup> are the coefficients. **Table 2** shows the estimated coefficients and the Standard Error (SE). It shows the t-statistic for each coefficient to test the null hypothesis of the coefficient being zero, given the other estimators in the model. The p-value of the F-statistic for the hypothesis of the coefficient being zero is also shown.

The statistics for the F-test on the regression model vs. constant model, showed significance of the model (F-statistics 15.1, p-value = 4.36 ∗ 10−<sup>9</sup> ). Standard deviation of the error distribution was 4.01 (weeks) and R-squared and adjusted Rsquared were 60 and 56%, respectively.

FHRV parameters were also used to estimate the GA. The following regression model for all healthy fetuses, without excluding the cases with low quality signals, were obtained using stepwise regression:

Estimated GA = b<sup>0</sup> + b1mRR + b2SDRR

where b0, b1, and b<sup>2</sup> are the coefficients. Therefore only the mean and standard deviation of fetal RR-intervals significantly contributed to the model. **Table 3** shows the estimated coefficients and SE of the coefficients. It shows the t-statistic for each coefficient to test the null hypothesis of the coefficient being zero, given the other estimators in the model. The p-value of the F-statistic for the hypothesis of the coefficient being zero is also shown.

The statistics for the F-test on the regression model vs. constant model, showed significance of the model (F-statistics 6.08, p-value = 0.004). However the standard deviation of the error distribution was 5.55 (weeks) which was larger than the SD

TABLE 2 | Results of Stepwise regression using valve intervals, including estimated coefficients (a<sup>0</sup> , a<sup>1</sup> , ..., a<sup>5</sup> ) and Standard Error (SE) of the coefficients, t-statistic and p-value for the F-statistic of the hypothesis of the coefficient being zero.


The model was obtained based on the parameters in milliseconds and GA in weeks. Other results include: F-statistics: 15.1, p-value = 4.36\*10−<sup>9</sup> , Standard deviation of the error distribution: 4.01 (weeks), R-squared: 60%, adjusted R-squared: 56%.

TABLE 3 | Results of Stepwise regression using FHRV parameters, including estimated coefficients (b<sup>0</sup> , b<sup>1</sup> and b<sup>2</sup> ) and Standard Error (SE) of the coefficients, t-statistic and p-value for the F-statistic of the hypothesis of the coefficient being zero.


The model was obtained based on the parameters in milliseconds and GA in weeks. Other results include: F-statistics: 6.08, p-value: 0.004, standard deviation of the error distribution: 5.55 (weeks), R-squared: 18%, adjusted R-squared 15%.

for the model with valve intervals, and R-squared and adjusted Rsquared were only 18 and 15%, respectively, which were smaller than those of the model with valve intervals.

Using leave-one-out cross-validation, the mean absolute difference between the CRL-estimated GA and the GA estimated from the proposed model was found to be 5.1 weeks using FHRV parameters and 3.8 weeks using valve timing intervals and 4.2 weeks when all parameters combined. When attempting to select a combined model, none of the FHR parameters were selected and therefore did not provide any additional value or increase the GA estimation accuracy. In the leave-one-out process, a new regression model is obtained when excluding each case. Therefore the model from the combined parameters were not necessarily the same as the model obtained based on all cases, thus the mean absolute errors were different for the methods using valve intervals and combined parameters.

# 3.2. The Effect of Signal Quality on GA Estimation

The absolute error which was calculated using leave-one-out approach, was not significantly correlated with the quality scores while controlling for GA and FHR:


However the error was decreased by applying a threshold on the quality of the signals and excluding the cases with low quality score. **Figure 3** shows the changes in the absolute error for various threshold values of fECG and DUS signal quality scores. Minimum absolute error 2.7 weeks was obtained when only the cases with fECG quality score > 0.4 and DUS quality score > 0.3

(22 cases in total) were considered. Applying higher thresholds would result in exclusion of more than 60% of cases. **Figure 4** shows the estimated GA using cardiac valve timings compared to the GA based on CRL as a gold standard for 22 fetuses with fECG quality score > 0.4 and DUS quality score > 0.3. This figure also shows the 1.96 × SD of the error as the 95% limits of agreement between the estimated age by cardiac valve intervals and CRL.

Furthermore, the absolute error for estimating the GA based on the FHR-related parameters could be reduced to the minimum of 4.7 weeks by applying the threshold of 0.32 on the fECG quality.

### 3.3. Changes of the Estimation Error with GA and FHR

The correlation of the error (estimated GA using valve intervals - CRL-based GA) using leave-one-out cross-validation, with GA as well as FHR was calculated and the results are summarized as follows:


• The error of GA estimation based on the valve timings was inversely correlated with FHR, when controlling for GA and quality scores (r = −0.325, p = 0.017). However, the absolute value of error was not significantly correlated to FHR (r = 0.1226, p = 0.3771).

## 3.4. Estimation of the GA for Abnormal Cases

The regression model based on the cardiac valve intervals of the healthy fetuses was used to estimate the GA of the fetuses with various abnormalities and arrhythmia. **Table 4** shows the GA estimated by regression and CRL for each case.

The table shows that arrhythmia (for cases 1–5) results in 5 weeks or less error in estimating the GA using heart valve intervals, while the FHR-based model failed to estimate the GA for bradycardia and arrhythmia case. **Figure 5** shows the estimated GA using cardiac valve timings vs. the GA based on CRL as a gold standard for the abnormal cases, compared to the 95% Confidence Interval (CI) for healthy cases (as shown in **Figure 4**). The specific abnormality types which resulted in estimated GA being outside the 95% CI are specified. The GA estimation using valve intervals clearly fails for some types of heart abnormalities such as ASD, VSD, SA and AV block (cases 19–22), due to their influence on opening and closing of the heart valves. FHR based model also failed for those anomalies as well as for the case with Premature Atrial Contraction (PAC) which was correctly estimated by the valve interval-based model.

## 4. DISCUSSION AND CONCLUSION

In this paper a new approach is proposed for estimation of the GA using fetal cardiac valve intervals. These intervals were estimated by a fully automated method from the raw recordings, therefore is less affected by human errors compared to sonography or LMP methods. Furthermore, the apparatus used to obtain the valve


TABLE 4 | Comparison of the estimated GA (weeks) using regression model based on the cardiac valve intervals with the gold standard (identified by CRL), for the cases with various types of anomalies and arrhythmias.

†Marks the estimated GA > 42 weeks, where the abnormal condition affects the valve intervals or FHR, therefore the regression model fails to estimate the GA correctly. The difference (estimated GA - CRL GA) for the cases marked with \* is outside the confidence interval of the error of GA estimation for healthy cases.

SSS, Sick sinus syndrome; WPW, Wolff-Parkinson-White syndrome; PAC, Premature Atrial Contraction; TOF, Tetralogy of Fallot; VSD, Ventricular Septal Defect; PA, Pulmonary Atresia; MS, Mitral Stenosis; ASD, Atrial Septal Defect; CDH, Congenital Diaphragmatic Hernia; CA, Chromosomal Aberration; AV block, Atrioventricular block; SA, Single Atrium; CAV, Cardiac Allograft Vasculopathy; CAVC, Common Atrioventricular Canal; PS, Polysplenia Syndrome; CHD, Congenital Heart Disease; NIHF, Nonimmune Hydrops Fetalis; TTTS, Twin-to-Twin Transfusion Syndrome Donor (Allan et al., 2000; Murray, 2007; Abuhamad and Chaoui, 2012; Allen et al., 2013).

intervals are easier to handle and require less skill to operate compared to standard sonography techniques. Although in this work we used non-invasive simultaneous recording of 1-D DUS signals and fECGs, the latter are only required for estimation of EDT which depends on the onset of QRS complex. It is also possible to use only 1-D DUS signals to obtain ICT, VET, IRT, and VFT using an automated technique that we proposed in previous work (Marzbanrad et al., 2013b). Since a 1-D DUS device can cost as little as \$17 and can be performed by nonexperts with limited training, it can be used to estimate the GA in resource limited settings (Stroux et al., 2014).

Although we have not directly evaluated the ANS and its relationship with valve intervals, the complex interplay between autonomic control of the heart and cardiac mechanics characterized by the valve intervals, has been previously reported in literature and is consistent with the results of this study (Berntson et al., 1994; Cacioppo et al., 1994a,b; Di Rienzo et al., 2013). According to the studies on adult cases, PEP is attributed to the sympathetic nervous system effect on the heart (Cacioppo et al., 1994a). As shown in **Figure 1**, the PEP, which is the duration from Q-wave to aorta opening, is comprised of two intervals: the EDT, which is the Q-wave to mitral closing interval, and the ICT, which is mitral closing to aorta opening interval. Our results show that not only do the ICT and EDT contribute to the GA estimate but their interaction is also a significant contributor. According to the results of stepwise regression in **Table 2**, VFT was also selected as a contributing term to the estimate the GA. Although less emphasis has been placed on

fetal VFT than other intervals, both in the literature and clinical practice, studies on adults found that VFT is controlled by both sympathetic and parasympathetic activity (Pinsky, 2005; Frazier et al., 2008; Khandoker et al., 2016). Therefore fetal development can be assessed by the ICT, the EDT and the VFT, as well as their interactions which evolve concomitantly with the changes in sympathetic and parasympathetic activities during fetal maturation. As discussed earlier, fetal autonomic brain age can be assessed using FHRV parameters (Van Leeuwen et al., 2003; Hoyer et al., 2013). The results of this current work demonstrated that a new method based on valve intervals outperforms the FHR-based method in estimating the GA although only time and frequency domain parameters and non identical populations were used. Furthermore, the valve interval method was less influenced by arrhythmias, particularly bradycardia, as shown in **Table 4**. FHR is also influenced by other factors such as behavioral states of the fetus and maternal physiological and psychological conditions, particularly in the second and third trimesters (Mantel et al., 1991; Monk et al., 2000; Ivanov et al., 2009; Marzbanrad et al., 2015b).While FHR might change according to those factors, this may not necessarily affect the estimation of GA by cardiac intervals. Among the cardiac valve intervals which were found contributing to estimation of gestational age, only VFT was correlated with FHR; the correlation of beat-by-beat fetal RR-intervals with ICT, EDT and VFT across 57 fetuses, were (−0.03 ± 0.13), (−0.02 ± 0.10) and (0.50 ± 0.21), respectively.

We note some limitations of the current study; first, the recordings used in this study are short and may not thoroughly represent the FHRV patterns which are used to evaluate the fetal functional brain development (Hoyer et al., 2013). Longer recordings would enable a better comparison of the effectiveness of valve intervals vs. FHRV patterns to assess the development of autonomic control. Further investigation using longer recordings is also recommended to be able to assess the influence of behavioral states and heart rate patterns on the valve intervals. We also acknowledge the recommendation of 5 min ECG for HRV analysis particularly for nonlinear measures. However, short term FHR variability e.g., the variation of beat-to-beat intervals for adjacent 3.75 s-epochs averaged over 1 min has been shown promising for monitoring of fetal development and surveillance of IUGR (Serra et al., 2008, 2009). On the other hand, the fact that most of the nonlinear HRV measures require at last 5 min of heart rate recording, further highlights an advantage of our proposed approach based on cardiac valve intervals, over the FHR-based approaches for GA estimation. Different from FHR-based approaches, the cardiac valve intervals require significantly less recording and measurement time to acquire a reliable estimation of the intervals to assess the fetal development. In clinical practice using ultrasound imaging, the parameters such as VET and PEP can be estimated by averaging over 30 s, and have shown to be well correlated with the gestational age (Mensah-Brown et al., 2010; Cruz-Martinez et al., 2012). This significantly reduces the examination time and the discomfort for the mothers.

The second limitation of our study is that our patient population did not include growth-restricted fetuses. To fully test our proposed technique, it would need to be evaluated on such a population. Although we have not evaluated the cardiac valve intervals for the fetuses with growth issues, this was studied in the literature, for young children (Alkon et al., 2003; van Deutekom et al., 2016). Van Deutekom et al., have shown that both birth weight and conditional height gain were independently associated with PEP, but not with Respiratory Sinus Arrhythmia (RSA).They discussed that as a shorter PEP indicates higher cardiac Sympathetic Nervous System (SNS) activity, this finding suggests that children with low birth weight have increased SNS activity compared with normal-birth-weight children. They associated the increased infant height gain with decreased SNS activity. Similar results specially for female children were observed by Feldt et al. (2011). Although these studies were not on fetuses, it is consistent with our result, showing that it might be extended to the fetal period. Furthermore, our study proposes a new physiological growth estimation method for healthy population, as the first step in identifying fetal development abnormalities. Furthermore, we studied fetal cardiac anomalies and arrhythmias which may confound the evaluation of GA, and obscure the potential detection of growth abnormalities. The third limitation is that the evaluation of the estimated GA using the proposed method for abnormal conditions, shows that the method fails to estimate the GA in presence of some, but not all, heart anomalies. This was particularly noted for the anomalies that affect operation of the valves. From another perspective, these results show that the valve intervals could be used to detect these anomalies that resulted in unrealistic GA estimates or large errors, as also investigated in our recent studies (Marzbanrad et al., 2014b; Khandoker et al., 2016). However, more cases with a large variety of anomalies are required for a more rigorous evaluation of their influence on the estimated GA.

Although in this study we have not considered the fetal gender, it might have an influence on the fetal growth and development. While some studies found no significant differences between male and female FHR during the first and second trimester (Neiger et al., 2004; McKenna et al., 2005), different intrapartum FHR patterns have been reported for two genders (Porter et al., 2016). Another study on term fetuses just before the labor, reported significantly lower values of most linear HRV measures for female fetuses compared to male fetuses, in both IUGR and control groups, as well as higher entropy indices in the control group (Gonçalves et al., 2013). Apart from FHR, the cardiac function was also previously investigated for male and female fetuses and no significant differences in cardiac function were found for different genders, except tricuspid valve E-wave velocity/time velocity integral for the entire tricuspid valve inflow (E/TVI) and pulmonary valve Acceleration Time (AT) (Clur et al., 2011). Based on the literature noted above, further investigations are required to study the gender-specific differences for GA estimation using valve intervals and FHRV. However, we note that in our study, the assumption was made that gender-determining technology would not be available. In other words, we aimed to make a system that could identify IUGR based upon the one dimensional Doppler only. In this context, the inclusion of gender into the algorithm would reduce the applicability of the approach we present.

Results of this study demonstrate that, for acceptable quality DUS and fECG recordings (determined automatically), the average error in GA estimation can be as low as 2.7 weeks, which is comparable to existing expert-driven methods. This proposed approach to GA estimation could be also improved with more accurate methods of quality assessment for 1-D DUS and fECG signals. It should be also noted that the error was obtained by comparing the estimated GA to the CRL estimates as gold standard, while a more accurate way would be to prospectively enroll the study subjects prior to conception and confirm the day of conception. The CRL is however subject to error (95% confidence interval of around 10 days), particularly in case of pathologies or unsuitable positioning (Grange et al., 2000; Callen, 2011). Furthermore the valve intervals have the advantage that they reflect physiological development of the fetus which is not completely aligned with the physical growth of the fetus. As discussed earlier, the sonography methods can be affected by genetic variations, such as the head shape, positioning of the fetus and pathologic conditions. While the error of the ultrasound-based GA predictors increases with gestational age (Caughey et al., 2008; Falatah et al., 2014; Al-Amin et al., 2015), according to our results, the error of the valve-based method does not change with gestational progression. The accuracy of the FHR-based method even increased with advancing gestation. Therefore our proposed physiological measures can be used in second and third trimesters, when the ultrasound imaging measures have a low accuracy and fail to detect abnormal growth, particularly in late gestation (Al-Amin et al., 2015). Overall, the proposed technique can be used as a measure of the physiological development and an adjunct in estimating the GA where ultrasound methods are unavailable or inadequate due to pathologies, unsuitable positioning, lack of skilled ultrasound operators, or other technical issues.

In conclusion, we proposed a novel and automated method for estimation of the GA, which could be performed using low cost, easy to operate devices that requires lower skills/training compared to sonography methods. In contrast to the sonography methods that are based on the physical growth, our proposed method provides assessment of the fetal physiological development. Compared to CRL-based GA estimates as gold standard, our method resulted in 2.7 weeks error for acceptable quality of recordings and also outperformed the GA estimation by FHRV parameters. The GA estimation based on valve intervals was affected by certain heart anomalies which influence the performance of the valves, but less affected by arrhythmias. Remaining errors in estimating the GA could be used as a marker to detect fetal abnormalities. Considering that the valve intervals reflect the autonomic control of the fetal heart, the new method provides automated assessment of the fetal ANS development that could be independent of the fetuses' locations on the growth curve (since our measure reflects neural development and not physical size). As a result the method proposed in this work might provide indications of growth-related issues, such as IUGR, early in pregnancy and potentially lead to early interventions.

#### AUTHOR CONTRIBUTIONS

FM, GC, and AK conceived and designed the study. YK performed the clinical experiments and data collection, while YK, AK, and MP carried out the fetal ECG experiments at Tohoku University Hospital. FM, GC, and AK contributed to the development of techniques and analysis tools. The data was analyzed by FM who implemented the signal processing and quality assessment, feature extraction, gestational age estimation, statistical analysis, comparisons and analysis of the results. All authors Participated in the discussion and interpretation of the results. FM, GC, and AK contributed to the writing of the manuscript and all authors reviewed and approved the final manuscript.

## ACKNOWLEDGMENTS

This study was supported by an Australian Research Council Linkage grant (LP100200184) between the University of Melbourne, Tohoku University and Atom Medical Corporation

#### REFERENCES


in Japan. A part of the work was carried out under the Coordination, Support and Training Program for Translational Research, Ministry of Education, Culture, Sports, Science and Technology in Japan. This research was also in part through the Len Stevens visiting scholarship to Emory University and Georgia Institute of Technology, and GC is funded by the Eunice Kennedy Shriver National Institute of Child Health and Human Development (NICHD) and the Fogarty International Center at the National Institute of Health (Grant# 1R21HD084114-01). The authors would also like to thank the clinical support service team at Tohoku University for providing the fetal Doppler and fECG data. The current study extends our previous conference paper (Marzbanrad et al., 2016) by adding fECG and 1D-DUS signal quality assessment which improved the results, analysis of the results for different quality and fetal heart rate, and analysis for 30 additional cases with abnormalities.


of fetal valve motion identification. Comput. Cardiol. Conf. 42, 365–368. doi: 10.1109/cic.2015.7408662


growth with childhood autonomic nervous system activity and its mediating effects on energy-balance-related behaviours the abcd study. Int. J. Epidemiol. 45, 1079–1090. doi: 10.1093/ije/dyw236


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Marzbanrad, Khandoker, Kimura, Palaniswami and Clifford. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Ultrasound Imaging of Mouse Fetal Intracranial Hemorrhage Due to Ischemia/Reperfusion

#### Kenichi Funamoto<sup>1</sup> \*, Takuya Ito<sup>2</sup> , Kiyoe Funamoto<sup>2</sup> , Clarissa L. Velayo<sup>3</sup> and Yoshitaka Kimura<sup>2</sup>

*<sup>1</sup> Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai, Japan, <sup>2</sup> Graduate School of Medicine, Tohoku University, Sendai, Japan, <sup>3</sup> College of Medicine, University of the Philippines, Manila, Philippines*

Despite vast improvement in perinatal care during the 30 years, the incidence rate of neonatal encephalopathy remains unchanged without any further Progress towards preventive strategies for the clinical impasse. Antenatal brain injury including fetal intracranial hemorrhage caused by ischemia/reperfusion is known as one of the primary triggers of neonatal injury. However, the mechanisms of antenatal brain injury are poorly understood unless better predictive models of the disease are developed. Here we show a mouse model for fetal intracranial hemorrhage *in vivo* developed to investigate the actual timing of hypoxia-ischemic events and their related mechanisms of injury. Intrauterine growth restriction mouse fetuses were exposed to ischemia/reperfusion cycles by occluding and opening the uterine and ovarian arteries in the mother. The presence and timing of fetal intracranial hemorrhage caused by the ischemia/reperfusion were measured with histological observation and ultrasound imaging. Protein-restricted diet increased the risk of fetal intracranial hemorrhage. The monitoring of fetal brains by ultrasound B-mode imaging clarified that cerebral hemorrhage in the fetal brain occurred after the second ischemic period. Three-dimensional ultrasound power Doppler imaging visualized the disappearance of main blood flows in the fetal brain. These indicate a breakdown of cerebrovascular autoregulation which causes the fetal intracranial hemorrhage. This study supports the fact that the ischemia/reperfusion triggers cerebral hemorrhage in the fetal brain. The present method enables us to noninvasively create the cerebral hemorrhage in a fetus without directly touching the body but with repeated occlusion and opening of the uterine and ovarian arteries in the mother.

Keywords: antenatal brain injury, fetal intracranial hemorrhage, ischemia/reperfusion, intrauterine growth restriction, cerebrovascular autoregulation, ultrasound imaging

### INTRODUCTION

Intrapartum asphyxia is a major cause of acute antenatal ischemic stroke and later irreversible neonatal encephalopathy (Cowan et al., 2003). Khatri et al. previously elucidated the pathophysiology for this type of brain injury as a double cascade of events involving ischemia and reperfusion leading to the disruption of the blood brain barrier and eventually, hemorrhagic cerebral transformation (Khatri et al., 2012). Known maternal antenatal risk factors for ischemia/reperfusion injury include primiparity, infertility, infection, pre-eclampsia, gestational

#### Edited by:

*Raimond L. Winslow, Johns Hopkins University, United States*

#### Reviewed by:

*Christoph Eugen Hagemeyer, Monash University, Australia Caterina Guiot, University of Turin, Italy*

\*Correspondence: *Kenichi Funamoto funamoto@fris.tohoku.ac.jp*

#### Specialty section:

*This article was submitted to Computational Physiology and Medicine, a section of the journal Frontiers in Physiology*

Received: *20 March 2017* Accepted: *09 May 2017* Published: *24 May 2017*

#### Citation:

*Funamoto K, Ito T, Funamoto K, Velayo CL and Kimura Y (2017) Ultrasound Imaging of Mouse Fetal Intracranial Hemorrhage Due to Ischemia/Reperfusion. Front. Physiol. 8:340. doi: 10.3389/fphys.2017.00340*

**92**

diabetes, smoking in pregnancy, and maternal nutrition (Nelson and Lynch, 2004; Wu et al., 2004; Lee et al., 2005; Darmency-Stamboul et al., 2012).

Our recent prenatal programming studies centered on maternal protein restriction illustrate the inherent susceptibility of fetuses and neonates to this type of neural compromise due to intrauterine growth restriction (Velayo et al., 2010, 2014; Dong et al., 2014, 2015). Molecular evidence show both prenatal and postnatal adaptive responses to protein-restricted diets rendering offspring vulnerable to further damage. These findings reinforce the widely-accepted statute that maternal nutrition in pregnancy is critical in early fetal brain development.

Currently, the exact mechanisms by which intrapartum hypoxic-ischemic events lead to acute antenatal ischemic stroke have not been fully elucidated due to the need for time-dependent observation. In this study, we adapted our intrauterine growth restriction mouse model for exposure to ischemia/reperfusion injury to evaluate the presence and timing of fetal intracranial hemorrhage. The ischemia/reperfusion cycles were generated for a mouse fetus by occluding and opening the uterine and ovarian arteries in the mother. Before and during the treatments, time-lapse ultrasound imaging of the fetal brain was performed. In addition, hemorrhage in the fetal brain was determined by histological observation. Our results indicate that protein-restricted diet increases the risk of fetal intracranial hemorrhage caused by ischemia/reperfusion, due to a breakdown in cerebrovascular autoregulation, showing the hemorrhage likely occurs after the second ischemic period.

### MATERIALS AND METHODS

#### Experimental Animals

This study was carried out in accordance with the recommendations of Regulations for Animal Experiments and Related Activities at Tohoku University, Center for Laboratory Animal Research, Tohoku University. The protocol was approved by the Center for Laboratory Animal Research, Tohoku University. Virgin female C57BL/6N mice (CLEA Japan, Inc., Tokyo, Japan) about 5 weeks old were maintained under controlled lighting (12-h light cycle) and temperature (24◦C). These were divided into two diet groups: Normal (N) and Low Protein (LP) by feeding them AIN-93G and modified AIN-93G (Oriental Yeast Co., Ltd., Tokyo, Japan), respectively. They were allowed free access to food and water ad libitum. The LP diet consisted of only 48% of the protein and 94% of the calorie contents of the N diet as suggested in previous studies (Reeves et al., 1993; Ito et al., 2011a,b). After a 2-week acclimatization period, the female mice were time mated and inspected for vaginal plugs the following morning. Plug-positive females were then transferred to single cages and fed ad libitum. Food consumption and body weight were recorded every day. There was a total of 25 pregnant mice per diet group.

#### Ischemia/Reperfusion Procedure

We modified Magal's method (Magal et al., 1988) as follows: on day 17.5 of gestation, about 1–2 days before birth, pregnant mice were anesthetized by subcutaneous injection with Ketamine (Ketalar intramuscular 500 mg, Daiichi Sankyo Propharma Co., Ltd., Tokyo, Japan; 10 mg/kg) and Xylazine (Selactar, Bayer Yakuhin, Ltd., Osaka, Japan; 5 mg/kg). Maintenance of anesthesia was achieved using inhalational 0.5% isoflurane at a 200 ml/min air flow rate (Forane, Abbott Japan Co., Ltd., Tokyo, Japan). Each pregnant mouse was placed on a 38◦C temperature controlled stage. Maternal abdomen was depilated with commercial hair removal cream (Veet, Reckitt Benckiser Group plc, Slough, England, UK), and this was followed by an abdominal midline incision. The maternal uterine horns containing 2–8 fetuses were exposed, and warmed with ultrasound gel (Parker Laboratories, Inc., Fairfield, NJ, USA) which also enabled ultrasound transmission described later. One fetus was chosen for ischemia/reperfusion treatment. Interval occlusion and opening of maternal uterine and ovarian arteries every 5 min induced ischemia and reperfusion (**Figure 1**). The ischemia/reperfusion sequence was repeated three times and the whole process lasted a total of 30 min. The number of fetuses treated in each diet group was n = 25.

#### Ultrasound Imaging

Using a Vevo 2100 high-frequency, high-resolution ultrasound device equipped with a linear array transducer MS-550D (VisualSonics, Inc., Toronto, Canada), ultrasound B-mode imaging of a cross section of a fetus brain was performed before and all throughout ischemia/reperfusion intervals (**Figure 1**). The central frequency of the ultrasound was 40 MHz resulting in the axial resolution of about 40µm. The field of view was 7 mm in width by 10 mm in depth, and three different focal depths, 4, 6, and 8 mm, were set, simultaneously. One hundred B-mode images were stored at 50 frame/s, wherein each captured cineloop duration was 2 s. Moreover, before and after the ischemia/reperfusion treatment, three-dimensional (3D) ultrasound power Doppler imaging was performed to visualize concomitant blood flow in the area. The scanning length for the 3D imaging was set as 4.98 mm at the interval of 0.102 mm.

#### Histology

To evaluate the presence or absence of intracranial hemorrhage after ischemia/reperfusion treatments, the fetuses were delivered and whole brain samples were collected. Brain tissues were supercooled in dry ice for 30 min and stored at −80◦C. These were subsequently mounted using an optimal cutting temperature (OCT) compound (Tissue-Tek 4583, Sakura Finetek Japan Co., Ltd., Tokyo, Japan) and cut on a cryostat (Leica Cryostat CM3050 S, Leica Microsystems GmbH, Wetzlar, Germany) to obtain cross sections of 8µm thickness. Images

of brain sections corresponding to ultrasound B-mode images were captured using Leica CTR 5000 and DM 5000B microscope system (Leica Microsystems GmbH, Wetzlar, Germany).

#### Image Analysis

Changes of blood flow in fetal brains were quantified based on image intensity of ultrasound B-mode images. First, in order to suppress noise in the B-mode images, 100 images obtained at each measurement time point (see **Figure 1**) were averaged, with adjusting for image shifting due to maternal breathing and heartbeat as well as fetal movement. A customized image analysis program calculated the cross-correlation function between the first image and each subsequent image in the measurement period. Each image was shifted so as to maximize the function, and the time-averaged B-mode image at each time point was then created by averaging image intensities at each pixel. Subsequently, the summation of image intensity, I, in a circular region of interest (ROI) with 15-pixel diameters of 0.35 mm was obtained with a separate image-processor (ImageJ, NIH, Bethesda, MD, USA). The ROIs were set at four different locations: cingulate cortex, area 2 (Cg2), basal forebrain (BF), and bilateral caudate putamen (Cpu). Finally, for each ROI, the summation of image intensities, I, was normalized with that before the treatment, I0, and the variation of the relative value of I/I<sup>0</sup> was evaluated. A total of five N fetal brains without hemorrhage and five LP fetal brains with hemorrhage, all previously confirmed by above-mentioned histology, were analyzed.

#### Statistical Analysis

The number of hemorrhage-positive fetuses on histology for each diet treatment group were counted, and the difference between the two groups were evaluated by Fisher's exact test. The differences in intensity of ultrasound B-mode images between the groups were analyzed using a t-test. Statistical significance for each analysis was set at p < 0.05.

# RESULTS

### Histological Findings

A snapshot of fetuses with or without intracranial hemorrhage after ischemia/reperfusion treatments is shown in **Figure 2A**. A dark region, indicated by a white arrow, in the head of the righthand side fetus indicates a hemorrhage. The number of fetuses being positive/negative for intracranial hemorrhage determined by histological observation in each group is summarized in **Figure 2B**. LP fetuses had a 5-fold higher frequency of intracranial hemorrhage than N fetuses; intracranial hemorrhage was noted in 16% and 76% of N and LP groups, respectively. Microscopic observation revealed that bleeding-prone areas were mostly located in the outer peripheral potion of the lateral ventricles, namely Cpu, as observed in the dotted circles in **Figure 2C**.

#### Ultrasonographic Findings

Sequential ultrasound B-mode images of the same cross section in **Figure 2C** were successfully obtained during the

ischemia/reperfusion treatments as shown in **Figure 3A**. This achievement enabled us to quantify changes in blood flow by intensity analysis of the sequential images described later with **Figures 3B–D**. At the Cpu beside the lateral ventricles where intracranial hemorrhage observed in the microscopic image, increases of the intensity were recognized in the sequential ultrasound B-mode image.

Dotted circles in the microscope image indicate the hemorrhage sites.

Blood flow in N and LP fetal brains before and after the ischemia/reperfusion treatments are demonstrated in 3D ultrasound power Doppler images of **Figure 4**. No significant change was observed in the N fetal brain before and after the treatment. On the other hand, several main blood flows disappeared in the LP fetus as indicated by triangles in **Figure 4**.

#### Intensity Analysis Findings

Changes in blood flow within the fetal brain were evaluated by intensity quantification of ultrasound B-mode images with each five N and LP mice. Sequential ultrasound B-mode images of five N fetal brains without hemorrhage and five LP fetal

ones with hemorrhage, confirmed by histological observation using microscopy, were analyzed. The variations of the relative intensities in ultrasound B-mode images in the four ROIs are shown in **Figures 3B–D**. Note that the results of Cpu were obtained by combining 10 data sets of bilateral Cpu in five fetuses. Image intensities in all the ROIs of N fetal brains without hemorrhage were slightly increased with fluctuations due to ischemia/reperfusion treatments. In contrast, the image intensities of LP fetal brains gradually decreased in the Cg2 and BF ROIs, but increased in Cpu ROI, showing a large standard deviation. The relative intensities in the Cg2 and BF ROIs of LP fetal brains were significantly lower than those in N fetal brains at some time points after the second ischemic period. No significant difference was detected between the relative intensities in the Cpu ROI between N and LP fetuses, but the fluctuation of relative intensity in the Cpu ROI during the ischemic/reperfusion treatments were opposite between the two groups.

#### DISCUSSION

The present study provides further evidence that ischemia/reperfusion events lead to the evolution of antenatal stroke. Moreover, it supports the assertion that the frequency of fetal cerebral hemorrhage increases in conditions of maternal-fetal undernutrition. It is known that fetal undernutrition causes diseases or disorders after birth or in adulthood (Woodall et al., 1996). The maternal/fetal low-protein condition also contributes to fetal cerebral hemorrhage, showing

that the frequency of fetal cerebral hemorrhage was statistically increased as shown in **Figure 2B**. Most bleeding-prone areas were the outer peripheral potion of the lateral ventricles of Cpu as displayed in **Figure 2C**. This tendency of fetal intracranial hemorrhage corresponds to the reference (Volpe, 2001), which reviewed an increased likelihood of periventricular leukomalacia in the presence of intraventricular hemorrhage.

Hypoxic condition was created by occluding the uterine and ovarian arteries and holding blood at the placenta. It is reported that brain cells can endure 15-min oxygen deficiency, and the present method of three times 5-min occlusion of the maternal arteries corresponds to the circumstance. Former studies utilized an animal model of rabbit for cerebral hemorrhage by repeatedly occluding and opening the carotid artery. In the animal model, blood supply to the brain was physically interrupted, and was impractical to apply to investigation of fetal intracranial hemorrhage. The present method enables us to noninvasively induce cerebral hemorrhage in a fetus through maternal arterial occlusion and directly observe the physiologic phenomenon of fetal brain sparing though real-time ultrasonography.

The variations of the relative intensity of ultrasound B-mode images represented different tendencies between the ROIs and nutrient conditions. In all ROIs, the blood flow in the N fetus retained or rather slightly increased during ischemia period. This implies that the N fetus has cerebrovascular autoregulation to maintain a constant blood supply to the fetal brain during the ischemia/reperfusion. In contrast, in both Cg2 and BF ROIs, cerebral blood flow in the LP fetuses decreased in the ischemia period and slightly returned in the following reperfusion period (**Figures 3B,C**). Whereas, the blood flow in Cpu of LP fetuses increased and became unstable with a relatively large deviation after the second ischemic period (**Figure 3D**). Interestingly, fluctuations of relative intensity in the Cpu ROI were opposite between N and LP fetuses. In 3D ultrasound power Doppler imaging of an LP fetal brain (**Figure 4**), decreased blood flow in a large tributary of the anterior cerebral artery was clearly recognized (Dorr et al., 2007). From all these observations, it may be concluded that the LP fetal brain fails to maintain blood flow during conditions of ischemia/reperfusion due to autoregulatory malfunction.

This study revealed "where" and "when" the intracranial hemorrhage occurs. However, "why" it occurs should be further investigated regarding mechanical weakness of a fetal intracranial vasculature and cell death. Decrease of blood flow in LP fetus observed in this study can be caused by lessened heart rate and decreased cardiac output. Further investigation would be performed by monitoring both maternal and fetal heart conditions with electrocardiography. Prevention method for fetal intracranial hemorrhage is also future work.

# CONCLUSION

The presence and timing of fetal intracranial hemorrhage caused by ischemia/reperfusion injury were evaluated using an intrauterine growth restriction mouse model. Histological analysis and ultrasound imaging were performed on mouse fetuses undergoing ischemia/reperfusion cycles which were generated by occluding and opening the uterine and ovarian arteries in the mother. Protein-restricted diets increase the risk of fetal intracranial hemorrhage due to a breakdown in cerebrovascular autoregulation.

# AUTHOR CONTRIBUTIONS

TI and YK designed the study; all authors substantially contributed to the experiments; TI and KiF took care of mice with specialized diet; TI performed surgical treatments, KeF and KiF measured fetal brains with ultrasound equipment with supervision of YK; TI performed histological observation; KeF analyzed ultrasound B-mode images; and KeF, TI, and CV wrote the manuscript with contribution of all authors.

# FUNDING

Part of the work was carried out under the Collaborative Research Project of the Institute of Fluid Science, Tohoku University. YK acknowledges funding from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) Translational Research Network Program.

# ACKNOWLEDGMENTS

The authors are grateful to Prof. Toshiyuki Hayase at the Institute of Fluid Science, Tohoku University for providing us with the ultrasound imaging system for the present study.

# REFERENCES


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Funamoto, Ito, Funamoto, Velayo and Kimura. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

digital media

of impactful research

article's readership