|Home | About | Journals | Submit | Contact Us | Français|
We propose a multi-stage approach using Wavelet and Hilbert transforms to identify uterine contraction bursts in magnetomyogram (MMG) signals measured using a 151 magnetic sensor array. In the first stage, we decompose the MMG signals by wavelet analysis into multilevel approximate and detail coefficients. In each level, the signals are reconstructed using the detail coefficients followed by the computation of the Hilbert transform. The Hilbert amplitude of the reconstructed signals from different frequency bands (0.1–1 Hz) is summed up over all the sensors to increase the signal-to-noise ratio. Using a novel clustering technique, affinity propagation, the contractile bursts are distinguished from the noise level. The method is applied on simulated MMG data, using a simple stochastic model to determine its robustness and to seven MMG datasets.
In the clinical evaluation of pregnant women, the analysis of uterine contractions could be useful in the prediction of labor, especially for patients at high risk of premature delivery.
The uterine contractions are a result of complex electrophysiological phenomena. Studies have shown that the uterine myometrical activity is low throughout pregnancy and increases significantly during term or preterm labor (Garfield et al 1998).
To date, several different techniques have been developed to measure uterine contractions. The most commonly used clinical approaches are the intrauterine pressure catheter (IUPC) and the tocography (TOCO) (Smyth 1957).
The IUPC is an invasive method that measures the intrauterine pressure via a catheter (Devedeux et al 1993). A change in pressure inside the uterus is reflected by displacement of the fluid in the catheter.
The TOCO is the most commonly used non-invasive technique to measure the mechanical deflections produced by the uterine contractions. The tocodynamometer, a strain-gauge based measurement device, records the deflection of the maternal abdomen during a uterine contraction. The strain on the tensometric transducer, which is strapped on the patient's abdomen, is proportional to the strength of the contraction. The TOCO technique is widely used by physicians as it is simple and almost risk free to the mother and the fetus. The major drawback of this instrument is the susceptibility to maternal motion artifacts.
While, IUPC is more reliable and accurate than TOCO, it is an invasive procedure that requires the rupture of the amniotic membranes, thus limiting its use to patients in whom delivery is necessary. Due to the poor predictive power of the TOCO and the invasive nature of IUPC, neither technique has been beneficial in the prediction of preterm labor or the diagnosis of true labor at term.
Currently, two methods are employed to record the physiological activity of the uterine contractions: (i) the electromyography (EMG), recorded by electrodes attached to the abdomen, and (ii) a newly established method, the magnetomyography (MMG), based on the recording of the magnetic fields that correspond to electrical fields. These techniques measure the electrical/magnetic activity on the surface of the maternal abdomen, which is a result of a sequence of bursts or groups of action potentials that are generated and propagated in the uterine muscle.
The uterine EMG, which reflects the original process of muscle fiber excitation, measures the action potentials of the myometrial cells by means of internal and abdominal surface electrodes. Garfield et al performed simultaneous recording of the EMG activity directly from the uterus and the abdominal surface of rats (Buhimschi and Garfield 1996). They were able to show conclusively that the signals recorded from the abdominal surface correspond to those generated in the uterus, suggesting that similar techniques can be used in humans. The recording of the EMG activity on the human uterus using abdominal electrodes was also reported (Garfield et al 1998).
The EMG has a high temporal resolution. However, because of the differences in the conductivities of tissue layers, the EMG signals are often filtered during their propagation to the surface of the abdomen.
A noninvasive technique that measures the magnetic fields associated with the action potentials is the uterine MMG. First, MMG recordings of spontaneous uterine activity were reported by Eswaran and colleagues (Eswaran et al 2002).
The MMG recordings have important properties: (i) independent on tissue conductivity (ii) detection of the signal outside the boundaries of the skin without making any electrical contact with the body and (ii) they are independent of references, which ensures that each sensor mainly records localized sources.
The automatic detection of uterine contractions has been attempted by different studies. Radhakhrishnan et al attempted to detect uterine contractions in EMG using higher order zero crossing analysis (Radhakrishnan et al 2000). The authors discuss the feasibility of this technique and conclude that the first-order zero crossing is most suitable to discriminate between contraction and noncontraction segments.
Based on the assumption that the uterus remains quiescent throughout most of the pregnancy but close to the time of delivery its electrical activity increases, a generalized synchronization index, as an indicator to track the spatial patterns of uterine myometrial activity, was proposed (Ramon et al 2005). The authors were able to show that the synchronization indices and their spatial distributions described uterine contractions.
In a recent study, La Rosa et al used a multiple change-point estimator along with the K-means clustering algorithm to detect uterine contractions in selected regions of a 151 SQUID sensor array (La Rosa et al 2008). They compute the relative energy variation (RMS) and the first-order zero crossing to discriminate between the time segments with and without contractions. However, this method lacks threshold detection (La Rosa et al 2008).
The frequency content of the electrical burst activity corresponding to the uterine contractions changes with time (Garfield et al 2005). Hence a simple spectral analysis may not be adequate to capture the uterine contractile activity. In the present work, we propose a combination of Hilbert-wavelet transforms (time frequency analysis) and affinity propagation (AP) clustering algorithm (Frey and Dueck 2007) to identify the uterine contractions. While wavelet transform decomposes the MMG signals into different frequency bands, the Hilbert transform quantifies the power in each frequency band. After this step, we apply the clustering technique on the Hilbert amplitude and create the binary decision signal which marks the contractile period.
To perform the time-frequency decomposition, one can also use the short-time Fourier transform (STFFT). This technique employs the spectral estimation procedure. As we plan to identify the contractions in small time windows (30 s duration), one has to use a fast Fourier transform (FFT) length of at least 10 s to get a frequency resolution of 0.1 Hz to capture the low-frequency contractile activity. However, there will not be ample windows to get a good spectral estimation. On the other hand, if one decides to use a longer FFT window, one will get a good estimation of the spectrum at the expense of losing the desired frequency range of 0.1 Hz. To quantify the low-frequency contractions, we need to estimate the spectrum with good frequency resolution. As described above, this can be achieved only by compromising for the reliability in the estimation. Hence, using this technique it may not be possible to correctly characterize the spectral property of the signals; therefore, we opted to use the wavelet approach to identify the contractions.
Six healthy pregnant women, one of preterm and five of term, in gestational age between 36 and 41 weeks were investigated (a total of seven data sets). The uterine activity was recorded with the 151 channel MEG system called SARA (SQUID Array for Reproductive Assessment). The sensors spaced approximately 3 cm apart and covering an area of 1300 cm2 are arranged in a concave array that covers the maternal abdomen from the pubic symphysis to the uterine fundus and laterally over a similar span (Eswaran et al 2002). The system is placed inside of a three-layer magnetically shielded room to reduce the interference of external magnetic fields with the biomagnetic fields generated by the human body. Each recording session was designed to last 20 min at a sampling rate of 250 Hz. Based on comfort of the pregnant women, the recording on the preterm subject MMGP201 36w0d lasted only 15 min and on one term subject MMGT201 39w0d we performed two shorter duration recordings with a break in between. The resulting two data sets were analyzed and presented separately. The study was approved by the Institutional Review Board and all subjects signed the informed consent prior to participation.
The duration of a uterine contraction varies between 40 and 60 s (Devedeux et al 1993) with MMG activity typically preceding the contraction. To effectively capture the low-frequency activity, it is essential to remove artifacts such as breathing and maternal/fetal magneto-cardiogram (MCG). To achieve this, we downsample the data to 32 Hz and band pass filter it between 0.1 to 1 Hz. To suppress the maternal breathing artifact, we use a notch filter between 0.25 and 0.35 Hz. For all filtering purposes, we use the Butterworth filter with zero-phase distortion.
The preprocessed signals are partitioned in 30 s disjoint inspection windows. A small portion at the end of the signal which is not integer multiple of the inspection window is discarded from the analysis. Data xi(t) from the ith sensor are convoluted with second-order Daubechies discrete wavelet db2 to produce approximate and detail coefficients as follows:
where Wx(u, 2j) contains the wavelet coefficients (approximate cj (u) and detail dj (u)) at time u at scale 2j in the jth level of decomposition. Detail coefficients dj (u) can also be independently computed as follows:
where ψ(t) and ϕ(t) represent the mother and father wavelets, respectively, and ψ(t) is represented as the linear combination of ϕ(t).
In broad sense, the approximate coefficients correspond to the low-pass coefficients and capture the trend in the data whereas the detail coefficients correspond to the high-pass coefficients and keep track of the fluctuations in the data. Most of the wavelets frequency decomposes the signals in a dyadic fashion. Thus, at the (j + 1)th level of decomposition, the wavelet coefficients become half in number compared to the jth level. This is equivalent to downsampling the data into half compared to the previous level. Based on this argument, if we consider sf as the sampling frequency in Hertz, at the jth level of decomposition, we will have frequencies ranging between sf/2j+1 and sf/2j (this is analogous to band pass filtering of the signal in the frequency bands of sf/2j+1 and sf/2j). In wavelet analysis, the signal in a desired frequency band can be obtained by performing inverse wavelet transform using only the detail coefficients that correspond to the frequency band. We use nine levels of decomposition which allows us to focus on the frequency of the uterine contraction (0.1–1 Hz).
In addition to time-frequency decomposition property, the db wavelets also effectively remove the (linear or nonlinear) trends in the data. The order of the wavelets determines the order of the trends (first, second, etc) up to which it can remain orthogonal. Because the MMG data are trend free, the inclusion of higher order terms is not necessary and the db2 approach is sufficient (Hussain et al 2009).
Usually, the variance of the wavelet coefficients is computed as an equivalent to the spectral power and is used to characterize the underlying dynamics. However, as higher decomposition levels are reached, fewer coefficients will be available and as a consequence the spectral power cannot be reliably computed. To overcome this impediment, we reconstruct the signal using the wavelet coefficient (at each level of decomposition) and compute the Hilbert transform of the reconstructed signals in each band. The reconstructed signal x(t) and its Hilbert transform h(t) can be represented as a complex analytic signal, x(t) + i · y(t), where . In chosen 30 s inspection window, we compute the mean value of the amplitude of the analytic signal as an equivalent to the power of the signal. Let yj,s(k) represent the power of the signal from the sensor s at the jth level of decomposition, where k is the average of the inspection time window. In order to improve the signal level, in the next step we sum yj,s(k) from all the sensors and from all the frequency bands and obtain y(k).
To discriminate between the contraction and non-contraction periods, we partition y(k) using the AP algorithm. The inputs are the similarities s(i, q) = −yi − yq2 between data points yi and yq, set to negative squared error (Euclidean distance). Real-valued messages are exchanged between data points of y(k) until the optimal number of exemplars (equivalent to the cluster centroid in the K-means algorithm) and corresponding clusters emerge. One of the main advantages of the AP clustering algorithm is that the a priori specification of the number of clusters to be detected is not required (Frey and Dueck 2007).
In order to separate the clusters representing the contractions from clusters of non-contractile regions, we compute the signal-to-noise ratio (SNR) as follows:
A certain percentage of the SNR is defined as a threshold (see section 2.1) and the clusters above this threshold are quantified as contractile patterns and the rest as non-contractile patterns. A binary signal is constructed by defining one in the instances corresponding to the exemplars that exceeded the threshold and zero at all other instances.
A proper selection of the threshold is a crucial part as it will determine the sensitivity/specificity of our algorithm, and it mainly depends on how we decide to base our algorithm. If we want the algorithm to operate highly sensitive, we will quantify more bursts at the expense of identifying erroneous burst patterns. On the other hand if we want the algorithm to operate highly specific, we identify only the prominent bursts at the expense of losing a few real bursts which are in the noise region. In this study, we base our algorithm to operate highly specific, and to estimate a threshold for our study we propose the following stochastic model based on the second-order autoregressive process (AR2):
where a1 = 1.81 and a2 = −0.9 are the parameters of the AR2 process, and η is the Gaussian white noise. To estimate the two starting (initial) parameters for the second-order autoregressive model (AR2), we fitted, for a particular data set (MMGT202 39w0d), the MMG data from each sensor using the second-order autoregressive process. This resulted in 148 different realizations for each, first and second, parameter describing the process. In the next step, we took the median of all the realizations and used the resulting coefficients in our stochastic model.
We also understand the robustness of our approach using the model. We would like to note that the stochastic model we propose here is just to mimic the amplitude of the MMG signals by tweaking the parameters of the model. This model is just sufficient to understand the limitation of our approach as to distinguish contraction burst from non-contraction activity. However, this model cannot explain the different frequency characteristic that may be present in the measured MMG signals.
We simulate s(t) for a duration of 30 min with a sample frequency of 32 Hz. Further, we consider the following two variables in our model: (i) θ, which determines the percentage of occurrence of the bursts in the data and (ii) ρ, which modulates the amplitude to the desired level. As done in MMG data analysis, we divide s(t) into disjoint 30 s windows. In order to determine the burst locations, we generate uniform distributed random number for each window. If the random number generated for a window is less than or equal to θ, we insert burst in that window. To model the contraction, we separate the data from the window (in which we decided to insert the burst) and band pass filter the same in the frequency band of 0.25–0.5 Hz using the fourth-order Butterworth filter with zero-phase distortion. In the next step, we insert the filtered data in place of the original data. This procedure is repeated individually for all windows that are chosen for the insertion of bursts and denote the modified signal as s′ (t). We modify the SNR of s′ (t) as follows:
where z(t) represents the Gaussian white noise and var(·) denotes the variance.
In the first step, analysis of the stochastic model is discussed, with the objective of threshold selection, and is followed by the presentation of the MMG results.
For the model described in (5) we vary the two parameters (i) θ—which determines the percentage of occurrence of the bursts in the data and (ii) ρ—which modulates the amplitude to the desired level. θ runs from 0.1 to 0.9 in steps of 0.1 (i.e., the inserted windows cover from 10 to 90% of the signal length in steps of 10%) and ρ runs from 1.2 to 9 in steps of 0.2 which yields SNRs similar to those observed in the real MMG data. We processed F(t) using the same approach as performed for MMG data using the Hilbert-wavelet approach. For each value of θ, we synthesized 40 different realizations of the data by varying ρ between the specified range. In the next step, for each realization a threshold value was determined so that only true contractile regions were picked. That is, different threshold values were tried until the minimum/none of the spurious events were picked. This resulted in 40 different intermediate threshold values. In addition to that we also computed the SNR in that particular realization as performed for the MMG data (see equation (3)). After this step, the mean of the thresholds, denoted with T and the mean of the SNRs, denoted with S, were computed for each ρ value. To this end, we modeled T and S as an exponential decay using equation (6), solved for the coefficients numerically using the Newton–Raphson approach and obtained the coefficients a = 1.179 and b = −0.273:
Figure 1 presents the case in which the first parameter θ was set to 0.1 (i.e., 10% of the data (six windows) containing bursts) and the second parameter ρ was at its minimum value, that is, a low SNR is produced (figure 1(a)). In figure 1(b) we present the data simulated from the model for θ = 0.1 and ρ 9.
The instances at which contractions are inserted are shown as thin dashed lines. Additionally, the marker computed based on the manually chosen threshold value that yielded the highest specificity is plotted over the synthetic signal.
From the simple model (see equation (6)), we clearly note that for SNR = 1, one has to select all the exemplars that are greater than 90% of SNR (0.9*SNR) as a burst pattern. However, this cannot reflect the real scenario. Hence, if the SNR ≤ 1 we set the algorithm not to mark any burst regions and consider this scenario as the limiting case of our method. For MMG signals, the SNR obtained using equation (3) is substituted in equation (6) to compute the threshold and identify the burst regions.
For a typical data set, MMGt205_40w1d, we present in figures 2(a)–(d) the mean of the Hilbert-wavelet transform (in each 30 s inspection window) summed across all sensors in four different frequency bands that range from 0.0625 to 1 Hz. It can be observed that, despite the notch filter, the 0.25–0.5 Hz band is the highest active band. However, MMG activity can also be seen in the neighboring bands suggesting, for this particular case, a spread of contractions through the entire 0.125–1 Hz range. In addition, figure 2(e) depicts the integrated signal y(k) from all sensors and the four bands. Red diamonds represent the instances that exceed the threshold and received the contractile pattern label.
Figure 3 illustrates the estimated contractions in the seven data sets. For each data set, the contraction marker is plotted over the preprocessed signal, obtained from one representative sensor with the highest SNR.
Here we apply the AP approach on the integrated Hilbert amplitude from all the sensors and from all the frequency bands and construct the contraction marker on this signals using the above procedure. Hence the contraction identified in this approach is characteristic of the highest active frequency band.
Accurate and objective methods to predict the nature of labor will be of great interest for obstetricians. From a cellular level, the uterine contractions are results of the transmission of action potentials in the myometrial cells of the uterine muscle. This electrical activity is transmitted to the different regions of the uterus through the coupling of the myometrial cells by means of gap junctions. It has been postulated by Garfield et al that the gap junction increases during the labor in various species (Garfield et al 1977). Further, spectral analysis of EMG shows an increase in power around 0.7–0.8 Hz for subjects in labor compared to the non-laboring subjects (Garfield et al 2005). A global synchronization of the electrical activity over the entire uterus has been related to successful progress of labor and delivery of the fetus. In the past, uterine contractile activity has been successfully measured by means of MMG (Eswaran et al 2002). However, to further understand the uterine dynamics and to identify the contractions (bursts) in the MMG signals, a more sophisticated method would be required. As uterine contraction is not a continuous process, to correctly understand the uterine dynamics, the spectral analysis has to be done exactly in the time window where the contraction has occurred. This is a tedious process as one has to manually score the contraction patterns prior to the application of the spectral analysis. Instead, a time-frequency approach to localize the contractile patterns is attempted in this work using Hilbert-wavelet transforms.
Based on this approach, the MMG power is classified into signal and noise components by using the affinity propagation clustering technique. However, the definition of noise level is ill defined in the MMG data. Hence, based on the SNR obtained by processing the MMG using the algorithm, we proposed a stochastic model to understand the performance of the approach. The data simulated using the model were processed in the same manner as performed for the real MMG data, and the SNR was computed. To this end, an appropriate threshold was chosen to identify the contractions incorporated in the model. In this procedure, we have chosen the threshold to operate highly specific to identify the true contractile events in the model. Based on this simulation, we observed that the threshold used to identify the contraction decays exponentially with the SNR. To quantify this behavior, we fitted an exponential model and identified the parameters. We then used this model to compute the threshold for the SNRs obtained for the MMG datasets. As shown in figure 3, our proposed method identified the beginning and the end of the contractile burst activity. This is the first step in the analysis process of making this technique a clinically viable tool. Once the events have been accurately identified, we can further process for the frequency content of the burst, the burst duration, number of bursts per minute and spatial-temporal spread of the burst across the sensor space. Using spectral analysis Garfield et al have shown that the power spectral density of uterine EMG bursts in patients during active labor peaked at 0.71 + 0.05 Hz as compared to non-laboring term (0.48 + 0.03 Hz) patients (Garfield et al 2005). Also, the power spectral density peak values were comparatively low for patients not in labor with respect to patients in active labor. Further the number of burst per minute and spread of the burst across the sensors will increase as the pregnant mother approaches active labor. All these factors are potential indicators for the clinical diagnosis of labor. The proposed method will provide an accurate temporal evaluation of the occurrence of the contractile activity. Along with the assessment of spatial characteristics such as phase synchronization between pairs of sensors as well as by considering different bio- and physiological parameters, our method could contribute to a more accurate diagnosis of preterm labor. In further studies, we plan to correlate these parameters to the outcome measures such as time to active labor from last MMG recording and cervical dilation.
In a recent study, Hill et al used recursive partitioning to identify gestational age specific and threshold values for infectious and endocrine biomarkers to predict preterm delivery (Hill et al 2008). The results indicate that, according to gestational age, two different biomarkers are the most accurate in predicting preterm delivery within 48 h. Therefore, to predict preterm labor with the current approach and thus to achieve a practical use in the clinical setting, future studies should concentrate on correlation analysis between MMG outcomes of contractile activity and bio-physiological measures.
In this work, we have successfully identified the bursts in the MMG signals using Hilbert-wavelet transforms. Future studies will focus on (i) the power in different frequency bands to provide a better understanding of the electro and magnetophysiology of uterine activity and (ii) predicting the nature of the labor. Further, we would plan to apply this approach in a longitudinal setting (multiple recordings from same subjects at different gestation age) to monitor the frequency change of the burst patterns throughout the entire gestational age. The threshold selection procedure will be based not only on the SNR, but also on physiological variables (e.g. duration of the contraction) to minimize the effect of false positives created by motion artifacts. To understand the spatio-temporal organization of the uterine dynamics, phase synchronization analysis will be pursued in the burst regions identified using the proposed method.
The study was funded by NIBIB/National Institutes of Health 1R01 EB007264-01A2.