Assessment of Features between Multichannel Electrohysterogram for Differentiation of Labors

Electrohysterogram (EHG) is a promising method for noninvasive monitoring of uterine electrical activity. The main purpose of this study was to characterize the multichannel EHG signals to distinguish between term delivery and preterm birth, as well as deliveries within and beyond 24 h. A total of 219 pregnant women were grouped in two ways: (1) term delivery (TD), threatened preterm labor (TPL) with the outcome of preterm birth (TPL_PB), and TPL with the outcome of term delivery (TPL_TD); (2) EHG recording time to delivery (TTD) ≤ 24 h and TTD > 24 h. Three bipolar EHG signals were analyzed for the 30 min recording. Six EHG features between multiple channels, including multivariate sample entropy, mutual information, correlation coefficient, coherence, direct partial Granger causality, and direct transfer entropy, were extracted to characterize the coupling and information flow between channels. Significant differences were found for these six features between TPL and TD, and between TTD ≤ 24 h and TTD > 24 h. No significant difference was found between TPL_PB and TPL_TD. The results indicated that EHG signals of TD were more regular and synchronized than TPL, and stronger coupling between multichannel EHG signals was exhibited as delivery approaches. In addition, EHG signals propagate downward for the majority of pregnant women regardless of different labors. In conclusion, the coupling and propagation features extracted from multichannel EHG signals could be used to differentiate term delivery and preterm birth and may predict delivery within and beyond 24 h.


Introduction
Preterm birth, defined as birth before 37 completed weeks of gestation, is a leading cause of neonatal morbidity and mortality and has long-term adverse consequences for health [1]. However, the problem lies not only in preterm birth itself but also in threatened preterm labor (TPL), which is defined as the presence of uterine contractions with no or limited evidence of cervical change between 20 and 37 weeks gestation and can be difficult to distinguish from active preterm birth [2]. TPL is the most common cause for the pregnant woman to seek institutional delivery care and it involves prolonged hospitalization, unnecessary medical interventions, the associated increase in expense, and aggravated anxiety for the pregnant woman and her family. In practice, less than half of the pregnant women with TPL will give birth prematurely. Accurate diagnosis of preterm birth is therefore clinically important. In addition, assessment of the time to delivery (TTD) is beneficial for both the healthcare system and the family. Based on TTD, medical resources could be allocated rationally and prepared in advance. The unnecessary hospital visit will Various labor prediction techniques and measurements have been proposed, such as cervical length, the Bishop score, fetal fibronectin, and tocodynamometer [3,4]. None of these techniques has been demonstrated to diagnose TPL and assess TTD objectively and precisely. Over the years, electrohysterogram (EHG) has been gaining popularity as a promising and powerful new tool for characterizing the parturition process mainly for the diagnosis of preterm birth [5], prediction of imminent delivery [6], and monitoring of uterine contractions [7,8]. EHG represents the spontaneous myometrial bioelectrical activity in the form of intermittent bursts of action potentials that trigger the mechanical contraction of the uterus [9]. EHG signals could be measured non-invasively with electrodes on the maternal abdominal surface. Throughout pregnancy, EHG signals change from scarce and poorly coordinated at the early stage to more intense and synchronized as delivery approaches [10,11].
Although a majority of studies recorded multichannel EHG signals, they usually selected a single channel for analysis [12,13]. Many features have been extracted from EHG signals to recognize preterm birth. Time, frequency, and time-frequency features [14] such as root mean square, median frequency, peak frequency, energy distribution, etc., have been used to characterize EHG signals. Previous studies found that uterine activities are nonlinear processes that change with time, and nonlinear signal processing techniques could thus provide additional information on physiological changes during pregnancy and close to delivery. Correlation dimension, sample entropy, and Lyapunov exponent [5,15] have been applied to describe the nonlinear interactions between billions of myometrium cells. Further, multichannel EHG records were investigated using multivariate sample entropy and fuzzy entropy [1,[16][17][18] to quantify the underlying dynamical structural complexity of multivariate physiological systems by taking into account cross-channel dependencies and to obtain information on the delivery onset.
Studies on the propagation velocity and propagation direction of EHG signals have significantly increased in the last decade [19,20]. However, no consistent results were reported. Mikkelsen et al. found that the uterine contractions expressed by EHG signals propagate both in the downward and upward direction [21]. Lange et al. reported no single preferred direction of propagation for the contraction bursts [22]. Xu found that uterine contractions are more likely to spread toward the center of the uterus [23]. A comparison of these results is not straightforward due to the different approaches and data sources. However, previous literature on the uterus unanimously reveals a special complexity of its electrical propagation properties; further studies are therefore necessary to clarify the mechanisms of the uterine activity for the prediction of delivery.
This paper aims to propose algorithms to characterize the coupling and propagation of multichannel EHG signals collected by our custom recording device, differentiate term delivery and preterm birth, and predict the delivery within and beyond 24 h.

Materials and Methods
The overall flow chart of the proposed method is shown in Figure 1, including data acquisition, signal preprocessing, feature extraction, and statistical analysis. Briefly, EHG signals were collected clinically, and signal preprocessing was carried out to reduce the interference, downsample EHG signals, and remove the outliers. Six features between multichannel EHG signals were derived and analyzed statistically. Sensors 2022, 21, x FOR PEER REVIEW 3 of 18

Data Acquisition
A total of 219 pregnant women between 33 and 41 weeks of gestation participated in this study. All of them were required to sign an informed consent. The study was approved by the Local Ethics Committee of Peking Union Medical College Hospital (ZS-1453) and was conducted following the Declaration of Helsinki (1989) of the World Medical Association.
Using our custom device, eight monopolar EHG signals were recorded simultaneously for approximately 30 min with a sampling rate of 250 Hz. Figure 2a shows the placement of the electrodes during the signal recording. As shown in Figure2b, the electrodes M1 to M4 were placed at the bottom of the uterus, 3-4 cm above the navel; M5 and M6 were placed 3-4 cm below the navel; M7 and M8 were placed at the cervix, 6-8 cm below the navel. M1 and M4 were placed symmetrically at 6-8 cm to the left and right of the midline, respectively; M2 and M3, M5 and M6, and M7 and M8 were placed symmetrically at 3-4 cm to the left and right of the midline; the reference electrode R and the ground electrode G were placed at the left and right iliac bones, respectively.
In this study, three bipolar EHG signals (B1, B2, and B3) were obtained from six monopolar recordings for the following analysis, since this configuration not only largely reduces the amount of interference present in the monopolar EHG recordings but also intuitively determines whether the uterine contractions propagate upward or downward.

Data Acquisition
A total of 219 pregnant women between 33 and 41 weeks of gestation participated in this study. All of them were required to sign an informed consent. The study was approved by the Local Ethics Committee of Peking Union Medical College Hospital (ZS-1453) and was conducted following the Declaration of Helsinki (1989) of the World Medical Association.
Using our custom device, eight monopolar EHG signals were recorded simultaneously for approximately 30 min with a sampling rate of 250 Hz. Figure 2a shows the placement of the electrodes during the signal recording. As shown in Figure 2b, the electrodes M1 to M4 were placed at the bottom of the uterus, 3-4 cm above the navel; M5 and M6 were placed 3-4 cm below the navel; M7 and M8 were placed at the cervix, 6-8 cm below the navel. M1 and M4 were placed symmetrically at 6-8 cm to the left and right of the midline, respectively; M2 and M3, M5 and M6, and M7 and M8 were placed symmetrically at 3-4 cm to the left and right of the midline; the reference electrode R and the ground electrode G were placed at the left and right iliac bones, respectively.

Data Acquisition
A total of 219 pregnant women between 33 and 41 weeks of gestation participated in this study. All of them were required to sign an informed consent. The study was approved by the Local Ethics Committee of Peking Union Medical College Hospital (ZS-1453) and was conducted following the Declaration of Helsinki (1989) of the World Medical Association.
Using our custom device, eight monopolar EHG signals were recorded simultaneously for approximately 30 min with a sampling rate of 250 Hz. Figure 2a shows the placement of the electrodes during the signal recording. As shown in Figure2b, the electrodes M1 to M4 were placed at the bottom of the uterus, 3-4 cm above the navel; M5 and M6 were placed 3-4 cm below the navel; M7 and M8 were placed at the cervix, 6-8 cm below the navel. M1 and M4 were placed symmetrically at 6-8 cm to the left and right of the midline, respectively; M2 and M3, M5 and M6, and M7 and M8 were placed symmetrically at 3-4 cm to the left and right of the midline; the reference electrode R and the ground electrode G were placed at the left and right iliac bones, respectively.
In this study, three bipolar EHG signals (B1, B2, and B3) were obtained from six monopolar recordings for the following analysis, since this configuration not only largely reduces the amount of interference present in the monopolar EHG recordings but also intuitively determines whether the uterine contractions propagate upward or downward.  In this study, three bipolar EHG signals (B1, B2, and B3) were obtained from six monopolar recordings for the following analysis, since this configuration not only largely reduces the amount of interference present in the monopolar EHG recordings but also intuitively determines whether the uterine contractions propagate upward or downward. A total of 219 pregnant women were grouped according to term delivery (TD), threatened preterm labor with the outcome of preterm birth (TPL_PB), threatened preterm labor with the outcome of term delivery (TPL_TD), and the EHG recording time to delivery (TTD). The number of pregnant women is shown in Table 1.

Signal Preprocessing
EHG signal energy is mainly concentrated in the range of 0.1 to 3-5 Hz [24]; in this work, we performed a fifth-order Butterworth band-pass filter between 0.34 and 1 Hz [16] to eliminate the main interferences caused by motion, respiration, and cardiac electrical signals. Then, EHG signals were downsampled at 25 Hz to reduce the computational cost of the data analysis [24]. The median absolute deviation (MAD) was computed in a 120 s window with a sample point overlap.
where X i denotes the amplitude of the ith point and median(X) denotes the median amplitude of all the points within a window. If |X i − median (X)|> 3MAD, the X i is regarded as an outlier, which was replaced by the linear interpolation of its nearby non-outliers.

Feature Extraction
We preferred to analyze the whole window analysis rather than the EHG-burst analysis. This decision was motivated by the lack of robust tools to automatically identify EHGbursts in a raw recording which usually requires the supervision of experts, especially in early gestational ages with the relatively low signal-to-noise ratio. By contrast, the whole window analysis has been proven to provide relevant information for predicting preterm birth. In this work, a whole recording analysis with a window length of 120 s and 50% overlap was performed to characterize the EHG signals, which is a trade-off between computational cost and information loss [25]. We then computed the median value of all the analyzed windows as a representative feature of each recording. In total, six EHG features between multiple channels were extracted to reflect their coupling and information flow.

Multivariate Sample Entropy
Multivariate sample entropy [18,26,27] can reveal correlations present in multichannel data and provide a robust relative complexity measure for multivariate data. The main advantage of this algorithm is the implementation of joint information embedded in a multivariate vector of an analyzed signal.
We first normalize the multivariate data with Z-score to cater for multichannel EHG signals obtained from different positions. Then, the multivariate embedded reconstruction was based on the composite delay vector.
We define multivariable time series {x k, i } N i=1 , k = 1, 2, . . . , p, where p denotes the number of variates or EHG channels, and N denotes the number of samples in each channel. The composite vector is constructed as follows: where M = [m 1 , m 2 , . . . , m p is the embedding vector and T = τ 1 , τ 2 , . . . , τ p is the time delay vector. In this study, The multivariate sample entropy was estimated using the following procedure: Compute the distance between any two composite delay vectors as the maximum of the following form: • For a given composite delay vector and an assumed similarity threshold r, count the number of instances P i where d[X m (i), X m (j)] ≤ r, i = j, then calculate the frequency of occurrence, and define a global quantity, • Extend the dimensionality of the multivariate delay composite vector from m to m + 1. For a given vector X m+1 (i), calculate the number of vectors Q i , such that d[X m+1 (i), X m+1 (j)] ≤ r, i = j, then calculate the frequency of occurrence, and define B m+1 (r) = 1 • The multivariate sample entropy is calculated as

Mutual Information
Mutual information is a measure of the amount of information that one random variable contributes to another variable. Formally, the mutual information is defined as follows: where X and Y are EHG signals from two channels, p(x, y) is the joint probability distribution of X and Y, and p(x) and p(y) are the marginal probability distributions of X and Y, respectively. MI is zero when X and Y are statistically independent; the larger the MI, the higher the correlation between X and Y.
where σ X and σ Y represent the standard deviation of X and Y, respectively, represents the covariance of X and Y. CC(X, Y) is the Pearson correlation coefficient of X and Y, ranging from −1 to 1. The greater the absolute value of CC, the stronger the correlation between X and Y. If CC = 0, it indicates that there is no linear correlation between X and Y. In our study, X and Y are two EHG signals.

Coherence
Coherence can reflect the linear correlation between two time series in the frequency domain. Briefly, coherence is an extension of Pearson's correlation coefficient in the frequency domain, and the coherence between signal X and signal Y can be obtained by normalizing the square of the cross-spectra by the auto-spectra using the following equation: where P xy (f) represents the cross power spectral of X and Y, and P xx (f) and P yy (f) are the power spectra of signals X and Y, respectively, for a given frequency f.

Partial Granger Causality
Partial Granger causality [28,29] takes into account partial correlation compared to traditional Granger causality to eliminate the exogenous links. We assumed the three-channel EHG signals X, Y, and Z, in which X and Y are represented by the joint autoregressive model: where X t is the value of signal X at moment t, X t − i and Y t − i are the values of signals X and Y at moment t − i, p is the order of the model, a i , c i , d i , and f i are the coefficients of the joint regression model, and ∈ 1t and ∈ 2t are the error of the joint prediction model at moment t. The variance/covariance matrix: To calculate the effect of Y t on X t in the context of the third time-series Z t , we extend the concept as follows: with variance/covariance matrix: Partial Granger causality is then computed as If PGC X→Y|Z > 0, the signal flows from X to Y; if PGC X→Y|Z < 0, the signal flows from Y to X. To determine the genuine direction of signal flow between X and Y, define the direct partial Granger causality.
If DPGC X→Y|Z > 0, the signal flows genuinely from X to Y, if DPGC X→Y|Z < 0, the signal flows from Y to X.
The direct partial Granger causality of the two EHG signals is as follows: If the majority of DPGC B1→B2|B3 , DPGC B2→B3|B1 , and DPGC B1→B3|B2 are greater than zero, we assumed the EHG signals propagate downward along with the uterus. If the majority are less than zero, then upward.

Transfer Entropy
Transfer entropy quantifies the amount of information transferred from one variable to the other. Importantly, transfer entropy is nonparametric and can capture nonlinear coupling effects. It has been shown that transfer entropy is a nonlinear extension of Granger causality.
Given two concurrently sampled time series X = {x 1 , x 2 , . . . , x N } and Y = {y 1 , y 2 , . . . , y N }, the transfer entropy from X to Y, termed TE X→Y , can be derived from conditional entropies as follows [30,31]: where p(y i , y is the conditional probability mass function of past observations of X and Y (the driving process) for Y (the target process), p(y i , y (l) i − t ) is the conditional probability mass function of past observations of Y (the driving process) for Y (the target process), i indicates a given point, τ and t are the time lags in X and Y, k and l are the block lengths of the past values in X and Y. In this study, τ = t = 2, k = l = 1.
Transfer entropy can exclude the common external environmental influence on twotime series. According to (24), if TE is infinitely close to 0, there is no obvious transfer relationship from X to Y. If TE > 0, information flows from X to Y. To determine the genuine direction of signal flow between X and Y, we defined the direct transfer entropy: Direct transfer entropy of the two EHG signals was calculated as follows: Among DTE B1→B2 , DTE B2→B3 , and DTE B1→B3 , if no less than two of them are greater than 0, the preferred direction of uterine contraction is downward. If no less than two of them are less than 0, the preferred direction is upward.
For both direct partial Granger causality and direct transfer entropy, we computed the percentage of pregnant women with the upward or downward EHG propagation in TD, TPL_PB, and TPL_TD, as well as for TTD ≤ 24 h and TTD > 24 h.

Statistical Method
Statistical analyses were conducted using SPSS 24 (IBM Corp., Armonk, NY, USA). The box diagram was used to describe the mutual information, correlation coefficient, coherence, direct partial Granger causality, and direct transfer entropy between any two channels. The six EHG features between channels were abnormally distributed according to the Shapiro-Wilk test and Q-Q plot. Therefore, we performed the Mann-Whitney U test to examine the features differences between the group of TD, TPL_PB, and TPL_TD, as well as the group of TTD ≤ 24 h and TTD > 24 h. In addition, a chi-square test was performed to examine the proportion of propagation direction. A significance level of p < 0.05 (two-tailed) was set for all analyses. Figure 3 shows the box diagrams of multivariate sample entropy, mutual information, correlation coefficient, and coherence, corresponding to TD, TPL_PB, and TPL_TD from left to right in each of the subplots. It can be seen from Figure 3a that multivariate sample entropy in both TPL_PB and TPL_TD is very significantly larger than that in TD (p < 0.01), indicating that the multichannel EHG signals from TPL are more irregular and more complex than TD. Figure 3b shows that mutual information in both TPL_PB and TPL_TD is significantly or very significantly smaller than that in TD (p < 0.05 or p < 0.01). Figure 3c shows that correlation coefficient of TPL_TD is significantly or very significantly smaller than TD (p < 0.05 or p < 0.01) between channels, and correlation coefficient of TPL_PB has similar results for B1&B3 and B2&B3. However, no significant difference was found in coherence of any two channels between groups (p > 0.05) see Figure 3d. In short, TPL_PB and TPL_TD have weaker inter-channel correlations than TD, and no significant difference was found between TPL_PB and TPL_TD (p > 0.05).

Statistical Method
Statistical analyses were conducted using SPSS 24 (IBM Corp., Armonk, NY, USA). The box diagram was used to describe the mutual information, correlation coefficient, coherence, direct partial Granger causality, and direct transfer entropy between any two channels. The six EHG features between channels were abnormally distributed according to the Shapiro-Wilk test and Q-Q plot. Therefore, we performed the Mann-Whitney U test to examine the features differences between the group of TD, TPL_PB, and TPL_TD, as well as the group of TTD ≤ 24 h and TTD > 24 h. In addition, a chi-square test was performed to examine the proportion of propagation direction. A significance level of p < 0.05 (two-tailed) was set for all analyses. Figure 3 shows the box diagrams of multivariate sample entropy, mutual information, correlation coefficient, and coherence, corresponding to TD, TPL_PB, and TPL_TD from left to right in each of the subplots. It can be seen from Figure 3a that multivariate sample entropy in both TPL_PB and TPL_TD is very significantly larger than that in TD (p < 0.01), indicating that the multichannel EHG signals from TPL are more irregular and more complex than TD. Figure 3b shows that mutual information in both TPL_PB and TPL_TD is significantly or very significantly smaller than that in TD (p < 0.05 or p < 0.01). Figure 3c shows that correlation coefficient of TPL_TD is significantly or very significantly smaller than TD (p < 0.05 or p < 0.01) between channels, and correlation coefficient of TPL_PB has similar results for B1&B3 and B2&B3. However, no significant difference was found in coherence of any two channels between groups (p > 0.05) see Figure 3d. In short, TPL_PB and TPL_TD have weaker inter-channel correlations than TD, and no significant difference was found between TPL_PB and TPL_TD (p > 0.05).

Comparison of Multivariate Sample Entropy, Mutual Information, Correlation Coefficient, and Coherence
(a)  Figure 4 shows the box diagrams of direct partial Granger causality and direct transfer entropy, corresponding to TD, TPL_PB, and TPL_TD from left to right in each of the subplots. Both direct partial Granger causality and direct transfer entropy of TPL_PB are significantly or very significantly smaller than TD (p < 0.05 or p < 0.01) for B1 → B2 and B2 → B3. Similarly, both direct partial Granger causality and direct transfer entropy of TPL_TD are very significantly smaller than TD (p < 0.01) between all channels. In short, the EHG signals flow less in TPL_PB and TPL_TD than TD, and no significant difference was found between TPL_PB and TPL_TD (p > 0.05).  Figure 4 shows the box diagrams of direct partial Granger causality and direct tra fer entropy, corresponding to TD, TPL_PB, and TPL_TD from left to right in each of subplots. Both direct partial Granger causality and direct transfer entropy of TPL_PB significantly or very significantly smaller than TD (p < 0.05 or p < 0.01) for B1 → B2 and → B3. Similarly, both direct partial Granger causality and direct transfer entropy TPL_TD are very significantly smaller than TD (p < 0.01) between all channels. In sh the EHG signals flow less in TPL_PB and TPL_TD than TD, and no significant differe was found between TPL_PB and TPL_TD (p > 0.05).    Table 2 summarizes the number and percentage of pregnant women with the upward or downward EHG propagation in TD, TPL_PB, and TPL_TD using direct partial Granger causality and direct transfer entropy, respectively. Regardless of TD, TPL_PB, and TPL_TD, EHG signals propagate downward for the majority of pregnant women in terms of both direct partial Granger causality and direct transfer entropy. No significant difference was found in the proportion of propagation direction between TD, TPL_PB, and TPL_TD with both direct partial Granger causality and direct transfer entropy (p > 0.05).  Figure 5 shows the box diagrams of multivariate sample entropy, mutual information, correlation coefficient, and Coh, corresponding to TTD ≤ 24 h and TTD > 24 h from left to right in each of the subplots. It can be seen from Figure 5a that multivariate sample entropy of TTD ≤ 24 h is very significantly less than TTD > 24 h (p < 0.01). Figure 5b shows that mutual information of TTD ≤ 24 h is significantly or very significantly larger than TTD > 24 h (p < 0.05 or p < 0.01). Figure 5c shows that correlation coefficient of TTD ≤ 24 h is significantly or very significantly larger than TTD > 24 h (p < 0.05 or p < 0.01). Figure 5d shows that coherence of TTD ≤ 24 h is significantly larger than TTD > 24 h (p < 0.05). In short, the closer to delivery, the lower the complexity of the EHG signal and the stronger the connection between channels.  Figure 6 shows the box diagrams of direct partial Granger causality and direct transfer entropy, corresponding to TTD ≤ 24 h and TTD > 24 h from left to right in each of the subplots. Both direct partial Granger causality and direct transfer entropy of TTD ≤ 24 h are significantly or very significantly larger than TTD > 24 h (p < 0.05 or p < 0.01), which indicates the closer to delivery, the stronger the information flow between EHG signals.  Figure 6 shows the box diagrams of direct partial Granger causality and direct transfer entropy, corresponding to TTD ≤ 24 h and TTD > 24 h from left to right in each of the subplots. Both direct partial Granger causality and direct transfer entropy of TTD ≤ 24 h are significantly or very significantly larger than TTD > 24 h (p < 0.05 or p < 0.01), which indicates the closer to delivery, the stronger the information flow between EHG signals. Table 3 summarizes the number and percentage of pregnant women with upward or downward EHG propagation in TTD ≤ 24 h and TTD > 24 h using direct partial Granger causality and direct transfer entropy, respectively. Regardless of TTD ≤ 24 h and TTD > 24 h, EHG signals propagate downward for the majority of pregnant women in terms of both direct partial Granger causality and direct transfer entropy. No significant difference was found in the proportion of propagation direction between TTD ≤ 24 h and TTD > 24 h with both direct partial Granger causality and direct transfer entropy (p > 0.05).  Figure 6 shows the box diagrams of direct partial Granger causality and direct transfer entropy, corresponding to TTD ≤ 24 h and TTD > 24 h from left to right in each of the subplots. Both direct partial Granger causality and direct transfer entropy of TTD ≤ 24 h are significantly or very significantly larger than TTD > 24 h (p < 0.05 or p < 0.01), which indicates the closer to delivery, the stronger the information flow between EHG signals.  Table 3 summarizes the number and percentage of pregnant women with upward or downward EHG propagation in TTD ≤ 24 h and TTD > 24 h using direct partial Granger causality and direct transfer entropy, respectively. Regardless of TTD ≤ 24 h and TTD > 24 h, EHG signals propagate downward for the majority of pregnant women in terms of both direct partial Granger causality and direct transfer entropy. No significant difference was found in the proportion of propagation direction between TTD ≤ 24 h and TTD > 24 h with

Discussion
In the present work, we obtained three bipolar EHG signals and performed a whole recording analysis with a 120 s sliding window and 50% overlap to characterize the EHG signals. Six EHG features were proposed to describe the coupling and information flow between multiple channels. Significant differences were found between TPL and TD for multivariate sample entropy, mutual information, correlation coefficient, coherence, direct partial Granger causality, and direct transfer entropy, and between TTD ≤ 24 h and TTD > 24 h.
Sample entropy was considered to be particularly appropriate for revealing EHG changes with pregnancy progression and delivery. Fele-zorz et al. found that univariate sample entropy significantly decreases as delivery approaches [32], suggesting the signal complexity decreases and its regularity increases. By contrast, in women with TPL under tocolytic therapy, no significant difference was found for univariate sample entropy to predict imminent delivery (TTD < 7/14 days vs. TTD ≥ 7/14 days) [6]. In this work, we assessed the structural complexity of multichannel EHG signal with multivariate sample entropy and showed significant differences between term delivery and preterm birth, and between TTD ≤ 24 h and TTD > 24 h. This result agrees with Ahmed who first characterized the interaction between the multivariate complex systems to successfully discriminate between women who finally delivered at term and those who did so prematurely [9]. However, we did not find any significant results between TPL_PB and TPL_TD when the recording was conducted far from delivery. The discrepancy between these two works may be because we conducted EHG recordings in women with TPL under tocolytic therapy. This latter has been shown to have a significant influence on uterine myoelectric activity [33], thus masking the subtle changes in uterine myoelectric activity through pregnancy. Previous studies showed the feasibility of predicting imminent delivery with a time horizon of 7 days in women with TPL. However, the performance of the model dropped dramatically if the time horizon was 14 days [6].
Mutual information and correlation coefficient significantly increased for TTD ≤ 24 h and TD group, suggesting stronger synchronization and association between multichannel EHG signals as delivery approaches. Similar results were also observed in the previous study [16]. This finding was physiologically related to the uterine myometrial cell excitability [15,34] and the formation of gap-junction as delivery approaches [35][36][37][38][39], which results in more intense and coordinated uterine electrical activity [11,39]. Again, we did not find any significant difference between TPL_PB and TPL_TD, which may be due to the tocolytic drug effect on the uterine myoelectric activity. In addition, we did not find any significant difference in coherence between TPL and TD, which may be because coherence was seriously impaired by instantaneous interactions. In this regard, the phase lag index [40] or weighted phase lag index [41] has been used to determine the true interactions in multichannel electroencephalography avoiding transient interactions, which may be more suitable for assessing the strength of cross-channel coupling.
In addition, we also assess the direction of EHG propagation which remains unclear in the literature. Some studies reported a predominantly downward propagation of the uterine electrical bursts in delivery women [42]. A study demonstrated that EHG bursts propagate both downward and upward, suggesting a multidirectional propagation pattern [21]. The multidirectional propagation was also reported by Escalona-Vargas using 151-channel magnetomyogram recordings [43], whereas other studies found no significant or preferred direction of propagation [22,44]. In this work, we used direct partial Granger causality and direct transfer entropy to comprehensively describe linear and nonlinear causality and information transfer at the organ level rather than at the cellular level or local uterine activity. That is, we used eight electrodes that were strategically positioned on the maternal abdomen covering practically the whole uterus. We found that the majority of EHG signals propagate downward to expulse the fetus, particularly, for the imminent delivery within 24 h. Our results agree with Garfield et al., who stated that uterine myoelectric activity presented a predominantly downward direction in women in the active phase of delivery using 3D vector myometrogram [20]. Using a 4 × 4 grid electrode with a relatively short inter-electrode distance, Diab proposed to use a nonlinear correlation coefficient and the index of general synchronization to determine the uterine myoelectric activity propagation and found that signals propagate in all directions but dominantly towards the cervix [45]. De Lau et al. also reported cases in which a downward propagating wave of uterine activity during a contraction was observed using a high-density grid of 64 electrodes. By analyzing the running cross-correlation of multichannel EHG records, Horoba found that the signal was generally delayed concerning that from the fundus for both physiological deliveries and threatened preterm labor [46]. Our results also agreed with Planque who found propagation in a descending direction in 87% of cases of women at the organ level [47]. We believe that the downward direction may occur in active phase delivery and at most a few days before delivery. For this reason, we only found dominant downward directionality in TPL_TD group and TTD ≤ 24 h, but not in TPL_PB and TD. Our results also stated that the propagation of uterine electrical activity does not show a preferential direction during pregnancy which seemed to be characterized by a highly unpredictable and potentially complex propagation pattern of individual spikes [48]. As far as we know, direct partial Granger causality and direct transfer entropy were first introduced in the present study to describe EHG propagation.
Due to the limitation of TPL sample size, more EHG signals from TPL will be collected in the next study to further validate the proposed algorithm. In this work, we found multiple multichannel features that could be used to determine delivery proximity. Further work is still needed to determine if this information is complementary or redundant to single-channel features for predicting term delivery and preterm birth. We analyzed the coupling strength and propagation direction in fast-wave high bandwidth which has been associated with cell excitability [49]. Future works could extend this analysis to fast-wave low bandwidth related to signal propagation to corroborate the downward propagation direction. We also noted that the features themselves do not have the prediction capability, and only the model trained by the features can be applied for prediction. However, the distinguishable features ensure the performance of the model. Our study attempted to explore the discriminable features which could be used to train classifiers in further study.

Conclusions
EHG is a very promising tool for monitoring uterine electrical activity with a wide range of applications. We extracted six EHG features between multiple channels from the bipolar recordings to distinguish term delivery from preterm birth, as well as deliveries within and beyond 24 h. Significant differences were found for these six EHG features between TPL and TD and between TTD ≤ 24 h and TTD > 24 h. We demonstrate that EHG multichannel features can distinguish different labors. Furthermore, stronger synchronization and association between multichannel EHG signals were exhibited as delivery approaches. Mostly, EHG signals propagated downward regardless of different labors.
In summary, the EHG features between multiple channels can provide coupling and propagation information to differentiate labors and facilitate the prediction of term delivery and preterm birth, and imminent delivery.