Shear Wave Dispersion as a Potential Biomarker for Cervical Remodeling During Pregnancy: Evidence From a Non-Human Primate Model

Shear wave dispersion (variation of phase velocity with frequency) occurs in tissues with layered and anisotropic microstructure and viscous components, such as the uterine cervix. This phenomenon, mostly overlooked in previous applications of cervical Shear Wave Elasticity Imaging (SWEI) for preterm birth risk assessment, is expected to change drastically during pregnancy due to cervical remodeling. Here we demonstrate the potential of SWEI-based descriptors of dispersion as potential biomarkers for cervical remodeling during pregnancy. First, we performed a simulation-based pre-selection of two SWEI-based dispersion descriptors: the ratio R of group velocities computed with particle-velocity and particle-displacement, and the slope S of the phase velocity vs. frequency. The pre-selection consisted of comparing the contrast-to-noise ratio (CNR) of dispersion descriptors in materials with different degrees of dispersion with respect to a low-dispersive medium. Shear waves induced in these media by SWEI were simulated with a finite-element model of Zener viscoelastic solids. The pre-selection also considered two denoising strategies to improve CNR: a low-pass filter with automatic frequency cutoff determination, and singular value decomposition of shear wave displacements. After pre-selection, the descriptor-denoising combination that produced the largest CNR was applied to SWEI cervix data from 18 pregnant Rhesus macaques acquired at weeks 10 (mid-pregnancy stage) and 23 (late pregnancy stage) of the 24.5-weeks full pregnancy. A maximum likelihood linear mixed-effects model (LME) was used to evaluate the dependence of the dispersion descriptor on pregnancy stage, maternal age, parity and other experimental factors. The pre-selection study showed that descriptor S combined with singular value decomposition produced a CNR 11.6 times larger than the other descriptor and denoising strategy combinations. In the Non-Human Primates (NHP) study, the LME model showed that descriptor S significantly decreased from mid to late pregnancy (−0.37 ± 0.07 m/s-kHz per week, p < 0.00001) with respect to the base value of 15.5 ± 1.9 m/s-kHz. This change was more significant than changes in other SWEI features such as the group velocity previously reported. Also, S varied significantly between the anterior and posterior portions of the cervix ( p = 0.02) and with maternal age ( p = 0.008). Given the potential of shear wave dispersion to track cervical remodeling, we will extend its application to ongoing longitudinal human studies.


INTRODUCTION
Preterm birth is the main cause of death in newborn children and its incidence rate is increasing [1,2]. A factor contributing to this problem is the lack of objective methods to describe the structural and functional changes that gestational tissues go through during pregnancy and labor [3].
One of the gestational tissues that suffers the most dramatic transformation is the cervix. The cervix is a cylindrical-shaped structure at the lower uterus that, during pregnancy, serves as a protective barrier between the developing fetus and the exterior. As pregnancy advances, the cervix softens and, close to delivery, shortens, ripens, and dilates to allow for vaginal parturition. In both term and preterm births, this transformation is achieved through a dramatic biochemical and structural transformation of the cellular and extracellular components of the cervical stroma [4].
The extracellular matrix of the cervix is composed of layers of semi-aligned collagen fibers immersed in a viscous matrix of long-chain molecules such as proteo-and glycosaminoglycans [5]. As pregnancy progresses, mature crosslinks holding together collagen fibrils are replaced by immature crosslinks. As a result, collagen organization is gradually lost and the cervix softens [4]. During the ripening stage, the concentration of sulfated glycosaminoglycans increases, particularly that of hyaluronic acid (HA) [6]. Due to the hydrophilic nature of HA, cervical hydration augments, further separating collagen fibrils [7][8][9][10][11][12][13][14]. All these changes alter drastically the mechanical properties of the cervix, such as its stiffness and viscosity [15].
During the last decade, various groups have investigated the use of Shear Wave Elasticity Imaging (SWEI) to objectively quantify changes in cervical stiffness during pregnancy and to assess the risk of preterm birth [16][17][18][19][20][21]. SWEI is an ultrasound-based method that generates shear waves in the tissue of interest using an acoustic radiation force impulse (ARFI) produced with a clinical transducer. The ARFI-induced shear wave motion is tracked with pulse-echo imaging techniques, allowing analysis of the propagation of the shear wave. If the medium is assumed to be linear, purely elastic, isotropic, and homogeneous, a simple relationship between the shear wave group speed V g and the tissue stiffness described by its shear modulus μ can be established: μ ρV 2 g , where ρ is the medium mass density [22]. Various studies have demonstrated a significant reduction of V g during pregnancy in ex vivo and in vivo studies in humans and animal models [17,20,[23][24][25]. However, these studies are limited by the aforementioned assumptions behind the link between the group velocity and the shear modulus, which ignore the effects that the complex extracellular matrix has on shear wave propagation.
The layered collagenous microstructure of the cervix and the presence of the highly viscous ground substance cause shear wave dispersion, i.e., the phase speed (the speed at which different frequency components of the shear wave move) varies with frequency. Shear wave dispersion has been demonstrated in other structurally complex tissues such as tendon and muscle [26,27]. Using a multi-scale simulation strategy with different viscoelastic models Peralta et al. [28] demonstrated significant shear wave dispersion in the cervix. More recently, our group observed that the slope of the phase velocity vs. frequency in ex vivo cervix samples from a non-human primate (NHP) model was 5.5 m/s-kHz (interquartile range: 1.5-12.0 m/s-kHz) [20]. In the same study, we also demonstrated that shear wave dispersion could be indirectly assessed by comparing the group speeds determined with particle velocity and particle displacement. The time derivative relationship between the particle velocity and the particle displacement makes the former more sensitive to higher frequencies, thus resulting in higher group speed values than those obtained with the particle displacement [29]. Interestingly, no differences in dispersion were observed between normal cervix samples and those in which ripening was induced with the prostaglandin Misoprostol. This was attributed to the large variance of shear wave displacements, probably related to the larger stiffness of the rhesus macaque cervix compared to the human one [20,30]. The quantification of descriptors of shear wave dispersion could provide quantitative imaging biomarkers that could be more sensitive than the group speed to changes in mechanical properties of the cervix.
Thus, in this study, we present evidence from a NHP model of the pregnant cervix that shear wave dispersion quantified with SWEI can be used as a quantitative imaging biomarker for cervical remodeling. We chose to work with the NHP model, in particular the Rhesus macaque, because of its similar maternalfetal physiology to humans. Furthermore, because Rhesus macaques are non-obligated quadripeds, the mechanical loading on the cervix is more similar to the human case [31][32][33][34][35]. Previous studies from our group have shown similar rates of shear wave reduction during pregnancy in humans and Rhesus macaques, thus showing the suitability of this model [21,24]. To overcome the limitations of high displacement variance in the Rhesus macaque cervix, we performed a pre-selection of two possible descriptors of dispersion, the ratio of group speeds from particle velocity/particle displacement and the slope of the phase velocity vs. frequency, combined with two noise reduction strategies. After the pre-selection stage, the descriptor with larger changes due to dispersion was applied to SWEI data from the mid-pregnancy and late-pregnancy NHP cervix.

SWEI-Based Descriptors of Dispersion
We investigated two possible SWEI-based descriptors of dispersion, one defined in the time domain, and another in Frontiers in Physics | www.frontiersin.org February 2021 | Volume 8 | Article 606664 2 the Fourier domain. The time-domain descriptor is the ratio R of group speeds computed with particle velocity v p over particle displacement u p : where v p can be determined from u p by: v p du p dt.
V g is usually estimated by tracking the time of flight of the maximum particle displacement or particle velocity away from the ARFI excitation [36]. In an approximately purely elastic material, both estimates of V g agree. However, this is not the case in viscoelastic media because the time-derivative in Eq. 2 acts as a high-pass filter in the frequency domain [29]. Thus, the more dispersive the medium is, the larger V g (v p ) is over V g (u p ) [20,37]. In other words, for a purely elastic material, R 1, and for a viscoelastic one, R > 1. This descriptor takes advantage of the various robust methods that have been proposed to quantify the group speed, such as the random sample consensus, or RANSAC, algorithm [38] or the Radon sum method [36,37] The Fourier-domain descriptor directly quantifies dispersion as the slope S of the phase velocity vs. frequency over the shear wave bandwidth [20]. To quantify S, we computed the particle velocity power spectrum V(k, ω), i.e., the square magnitude of the 2D-Fourier transform for the particle velocity v p (x, t), where k is the wave number and ω is the angular frequency [39]. To parameterize the frequency dependence of the phase speed V p (ω) we used a linear model defined by an intercept V p,0 and a slope dV p /dω. To find the values of the parameters of the linear fit, a Radon-sum strategy proposed by [40] was used. Briefly, curved trajectories are defined on the velocity power spectrum as function of the linear fit parameters. The values of the parameters are selected from the trajectory along which the sum of the spectrum values is maximized. The sum was done over the -6 dB spatial and temporal frequency bandwidth. Finally, we used S dV p /dω as the Fourier-domain descriptor of dispersion.

Pre-selection Stage
The pre-selection stage consisted of selecting the descriptor with the largest changes in dispersion in viscoelastic media over the estimation variance. To have control over the degree of dispersion, we used computational phantoms based on finiteelement models of Zener viscoelastic solids implemented on LS-Dyna (Version 3.2 β, Livermore, Software Technology Corp., Livermore, CA) following the method described by [41]. Dispersion was controlled by varying the viscous component of the Zener model (one spring in parallel with a series of a spring and a dashpot). Shear waves were created through fields of axial stresses with spatial Gaussian symmetry that simulated acoustic radiation force impulse stimuli produced by clinical ultrasound transducers [29]. The width of the Gaussian profile was 0.1 mm, based on a previously reported computational optimization [42]. The stress fields were applied during 0.1 ms. The center of the stress field was at a depth of 30 mm, corresponding the axial center of the simulated computational phantom [41].
The stress relaxation response of the Zener model has three parameters, the instantaneous shear modulus G 0 , the equilibrium shear modulus G ∞ , and the time constant b. The values of G 0 and G ∞ were fixed at 20 kPa and 2 kPa, respectively, and the value of b was varied from 3,000 s to 30,000 s, corresponding to changes in the viscous component from 0.6 to 6.0 Pa s. The temporal and spatial sampling rates of the simulated mechanical response (shear wave displacements) were 0.1 μs and 0.25 mm, respectively, for a total simulated time of 10 ms. For the shear wave frequencies that were simulated (peak temporal frequencies between 100 and 200 Hz as shown in Figure 1), there were at least 10 nodes per wavelength. To mimic the pulse repetitions frequencies used in commercial scanners such as the one used in the animal model study, the displacements were down sampled to 0.1 ms. The axial and lateral length was 60 mm and 20 mm respectively; additionally a 5-element thick perfect matching layer was added to all outward faces of the mesh to suppress reflections back into the analysis ROI. The analysis of shear wave propagation was done over a 5 mm axial extent centered at the ARF focus (30 mm deep), similarly to the analysis done in the animal model. Zero-mean Gaussian noise was added to particle displacements to analyze the compromise between changes of the descriptors due to dispersion and estimation variance [43]. The level of noise was controlled by a factor k that defined the standard deviation as a factor of the maximum displacement induced by the stress field in each simulation. The value of k was varied from 0 to 0.1.
The compromise between changes of each parameter and the estimation variance was quantified in terms of the contrast-tonoise ratio of each descriptor with respect to its value for the least dispersive medium. The CNR was mathematically defined as: where P η i is the value of either S or R for the material width the ith value of viscosity η, and P η min is the value of the same parameter for the least viscous medium. σ P η i and σ P η min are the standard deviations of these parameters. The value of the CNR for each parameter was computed for different combinations of viscosity and noise. Twenty independent realizations of the noise were added to the simulated shear wave fields to estimate the mean and standard deviation of the CNR.

Noise Reduction Strategies
Preliminary estimates of shear wave dispersion in the ex vivo NHP cervix [20] showed large variance possibly related to the high stiffness of cervical tissue compared to other soft tissues such as liver. To overcome this limitation and increase the sensitivity of dispersion descriptors, we implemented two noise reduction strategies.

Automated Low-Pass Filter
Because random noise dominates at high temporal frequencies [44], low-pass filters have been commonly used in previous applications of shear wave elasticity imaging on cervical tissue [24]. We implemented this strategy and designed an automated criterion to select the cutoff frequency. This allowed us to automatically adapt the filter to the noise level in each simulated condition. To this end, we first obtained the maximum projection of the 2D spectrum of the particle velocity to obtain the velocity power spectrum vs. temporal frequency. Then, the cut-off frequency was defined by the intersection of two lines fitted to the velocity power spectrum vs. temporal frequency within a bandwidth that contains the shear wave signal (100-500 Hz) and a high-frequency bandwidth that contains noise (1)(2)(3)(4)(5). This is graphically represented in Figure 2.

Singular Value Decomposition Filter
The second strategy consisted of a filter based on the singular value decomposition (SVD) of shear wave displacements. This technique has been used to reduce clutter in other ultrasound techniques such as Doppler imaging and functional ultrasound [45]. The SVD filtering strategy is based on the fact that noise components of the displacement field have low spatial and frequency correlation compared to the components caused by shear wave propagation. Thus, unlike the low-pass filter strategy defined above, the SVD filter considers both temporal and spatial information. Also, in contrast to previous uses of the SVD, we aim at maintaining highly correlated components of the displacement field. To apply this filter, the 2D displacements are sorted into a Casorati Matrix S with the spatial information in the FIGURE 2 | Projected particle velocity spectrum of a viscoelastic solid with 6.0 Pa s viscosity and a k 0.1 noise level (orange curve). The green and yellow lines are linear fits to the noisy spectrum at low (0.1kHz-0.5 kHz) and high frequencies (1 kHz-5 kHz), respectively. The intersection of both fits allows determination of the cutoff frequency f c . The blue curve shows the noise-free case to demonstrate the agreement of the fit over noisy data to the noiseless case. columns and the temporal information in the rows [45]. The Casorati matrix is decomposed into the product of three matrices, S UΔV p , where matrices U and V contain the eigenvectors of the covariance matrices SS p and S p S, respectively, and Δ is a diagonal matrix containing the singular values (square root of the eigenvalues of the covariance matrices) in decreasing order [46]. The first diagonal elements represent the contribution of the shear wave to the covariance (high contributions), while the last ones represent the contribution of random noise (low contributions). To filter noise, a displacement matrix S ′ is reconstructed removing the eigenvalues smaller than the cutoff eigenvalue, corresponding to the seventh eigenvector. The cutoff index value was defined by minimizing the tradeoff between bias and variance in the estimation of dispersion based on simulations. This tradeoff was quantified in terms of the Contrast-to-Total-Error-Ratio (CTER), defined as the ratio of the contrast of parameter P with respect to its value in the least dispersive medium over the combined root mean squared error RMSE: where RMSE was computed as: where P η,e is the expected value of the parameter for each degree of viscosity. The CNR of R and S for media with different dispersion (from different degrees of viscosity) and different noise levels where compared with and without each of the noise reduction strategies to identify the combination that maximized the CNR.

In vivo Application to Pregnant Rhesus Macaques
After identifying the combination of the dispersion descriptor and noise reduction strategy that maximized CNR, the combination was used to evaluate changes in the dispersion properties of the NHP cervix during pregnancy. To this end we used data from a longitudinal SWEI study on 18 pregnant Rhesus macaques acquired at weeks 10 and 23 of gestation, corresponding to mid and late pregnancy [21]. The study was performed at the Wisconsin National Primate Research Center and the protocol was approved by the Institutional Animal Care and Use Committee of the University of Wisconsin-Madison. Characteristics of the subjects are described in Table 1.
SWEI data was acquired from the NHP cervix using a transrectal scanning approach on a longitudinal view of the cervix. Data was acquired with a prototype catheter transducer on a Siemens Acuson S2000 ultrasound scanner (Siemens Healthcare, Ultrasound Business Unit, Mountain View, CA, USA). Technical details about SWEI implementation can be found in [21]. Shear waves were induced in two propagation directions (from the uterine end to the vaginal end of the cervix, and vice versa) and were tracked over a 6 × 5 mm 2 region of interest centered first on the posterior and then on the anterior portions of the cervix.
Shear-wave induced displacements were estimated from raw IQ data obtained from the system's Axius Direct Ultrasound Research Interface [47] with Loupas' correlation method [48]. Displacement estimates were rejected based on values of the correlation coefficient ( < 0.98) and filtered to eliminate low frequency biological motion using a second-order polynomial motion filter [40]. The first lateral millimeter of the displacement field was discarded due to ARF stimulus effects. In vivo data quality was assessed through three criteria ( Figure 3): 1. Evidence of a planar wavefront: Maps of the laterally-integrated particle displacement vs axial position and time must show no variation over the axial position (planar wavefront), and should gradually decrease over time mainly due to attenuation ( Figure 3A). Data not showing these characteristics were rejected ( Figure 3D). 2. Statistically significant particle displacements. Following a criterion proposed by [20], particle displacements were compared to those estimated from reference tracking frames acquired before the application of the ARF stimulus ( Figure 3B). Data was rejected if time-averaged displacements (pink curve in Figures 3B,E) were not larger than a threshold (red dashed line) defined as three times the standard deviation of the reference particle displacements (black dashed line) over more than half of the axial span of the ROI. For the rest of the data, we only processed displacements that were continuously above the threshold (right side of vertical green line in Figure 3E). 3. Data with high noise after filtering: Finally, we rejected all the filtered data with noise levels at high temporal frequencies above -6 dB of the maximum signal ( Figures 3C,F).

Statistical Analysis
We used a linear mixed effects model to evaluate the significance of the relationship between the optimized dispersion parameter, gestational age, and other experimental variables [49]. According to this model, the optimized parameter P can be written as a linear combination of each variable e i and of their interactions e i pe j : where P 0 is the intercept of the model, C i is the contribution of each variable to the model, and C i,j is the contribution of the interaction of each pair of variables. Table 2 defines the variables e i . We used a maximum likelihood strategy to fit the model using the Matlab's function fitlme. The model was applied in a sequential manner, discarding non-significant variables in each iteration. We report the final model that included significant variables only. A Shapiro-Wilk test was applied to the values of the optimized parameter grouped by gestation week and cervix portion in order to verify the assumption of normality.

Simulations
The CNR of both dispersion descriptors R and S increased with viscosity. In addition, we computed the values of the dispersion descriptors R and S (see description below) for a purely elastic medium (no viscosity), obtaining R 0.998 and S 0.01 (close to the values expected under no dispersion, i.e., R 1 and S 0). This demonstrates the absence of surface effects in shear wave propagation. Also, when no noise-reduction strategies were applied, R was found to be more robust to high levels of noise than parameter S. Figures 4A,B show the CNR of R and S, respectively, as a function of viscosity. Each curve corresponds to a different noise level. For low noise levels (k ≤ 0.02), the CNR of S  was higher than the CNR of R by a factor of six for the low viscosity material and a factor of 5.5 for the high viscosity material. However, beyond k 0.05 the CNR of S was negligible compared to the CNR of R.
For the low-pass noise filtering method, the automaticallyselected cut-off frequency increased as a function of viscosity and decreased as a function of the noise level. This is shown in Figure 5, where the cutoff frequency is plotted as a function of the noise factor k. Each curve corresponds to a different degree of viscosity.
For the SVD noise filtering method, the optimum singular value cutoff was chosen based on the maximum value of the CTER.   shows the CTER as a function of the singular value cutoff index for the material with 0.8 Pa s viscosity and the material with 6 Pa s viscosity. The CTER had a maximum value around the index seven for both materials. Thus, we selected a value of seven for the cutoff index. Figure 7 shows parameter S obtained after applying the lowpass filter (a) and the SVD filter (b), as a function of the noise parameter k and for various degrees of viscosity. With the low pass filter, changes in S as viscosity increases were not significant as indicated by the overlapping error bars. Furthermore, the highest value of S does not occur with the higher values of viscosity, suggesting ambiguity in the values of S that result from the use of this filtering strategy. This did not occur with the SVD filter, for which S increased monotonically with viscosity and changes with viscosity did not overlap. Because of this, the low pass filter was rejected.
Frontiers in Physics | www.frontiersin.org February 2021 | Volume 8 | Article 606664 9 After applying the SVD filter, the S parameter showed the highest CNR for each degree of noise and viscosity. Figure 8 shows the CNR of both R and S as a function of the noise factor without any filter and with the SVD filter. The combination of descriptor S and the SVD filtering for the highest level of noise resulted in 11.6 times higher CNR than R descriptor with the same filtering strategy. We applied the SVD filter to the estimation of descriptor S to data from the in vivo NHP study, zeroing singular values at and beyond index 7. Table 3 summarizes the mean and standard deviations of the S parameter at weeks 10 and 23 for the anterior and posterior regions of the 18 macaques. Empty cells correspond to cases in which no measurements passed the quality criteria described above. No standard deviations (ND) are reported for those cases in which only one measurement passed the quality criteria. Mean cervical thickens for each animal is also reported in Table 3. Figure 9 shows an example of the particle displacements (a,d), particle velocities (b,e) and dispersion curves (c and f) for data from different subjects at week 10 (top row) and week 23 (bottom row). In this example, the maximum values of the particle velocity propagate at a faster rate at early pregnancy than at late pregnancy, which indicates overall stiffer tissue in the former case. Dispersion curves show higher phase velocities vs. frequency at mid pregnancy than at late pregnancy.

In vivo NHP Application
In general, S decreased between weeks 10 and 23 in the anterior cervix portion (Wilcoxon p 0.03). Shear wave dispersion also decreased in the posterior portion, but the difference was not statistically significant (Wilcoxon p 0.051). Figure 10 shows box plots that group estimates of S for all the subjects that passed the quality criteria at weeks 10 and 23 for (a) the anterior and (b) the posterior portions of the cervix.
With exception of S values in the posterior side at week 23 (p 0.065), all data complied with the assumption of normality. Table 4 summarizes the results of the final LME model after sequentially rejecting non-significant effects. The most significant effect was that of the week of gestation (p < 0.00001), which resulted in a reduction of 0.37 ± 0.07 m/s-kHz per week from the initial value of 15.5 ± 1.9 m/s-KHz. Also, dispersion decreased with the mother's age (p 0.008) and was higher in the anterior portion (p 0.0009). Furthermore, significant correlations were found between the week of gestation and the cervical portion, and the maternal age and the cervical portion.

DISCUSSION
The goal of this study was to evaluate the changes in shear wave dispersion quantified through SWEI in an in vivo NHP model of the cervix between mid and late gestation. To this end, we performed a pre-selection of two SWEI-based descriptors of shear wave dispersion, R and S and noise reduction strategies using finite-element simulations of viscoelastic media. The most important findings of this study are: • The S parameter with the SVD noise reduction strategy showed the highest dispersion-based CNR compared to other combinations of parameter and noise reduction strategies. • The combination of the S parameter and the SVD filter led to observe a significant reduction of shear wave dispersion in NHP cervix between mid and late pregnancy. • As result of the application of the LME model, we found other significant effects that influence the value of the S parameter and included the mother's age and the cervix portion (anterior vs. posterior).
Based on a finite-element model, we observed that the frequency-domain parameter S in combination with the SVD filtering strategy produced the largest CNR among materials with different degrees of dispersion. The advantage of the SVD method over the automated low pass filter is based on its capacity to reduce both spatial and temporal contributions to noise. A disadvantage of this technique is the need to define an optimum cutoff value. A particularly complicated situation is how to treat frequency content due to bulk motion. In order to reduce the effects of such motion in the implementation of the optimized SVD strategy, we applied a second-order motion filter, as done by [40]. In the future, we will investigate the optimization of the selection of the singular value cutoff index to also reduce the contribution of bulk motion.
Other studies have applied SVD filtering in ultrasound studies like blood flow imaging [50]. However, the optimization of the cutoff value of SVD filtering for SWEI has not been widely studied. A particular case of SVD filtering in shear wave analysis can be found in geophysics applications. Chevrot and Girardin [51] used SVD filtering to reduce the noise in shear waves from seismic signals. In that study, a cutoff value was selected when a regular decrease of the singular spectrum was observed, which led to a cutoff index of 6. Here, we proposed a quantitative selection criterion considering the bias and variance in the estimation of the dispersion slope based on simulations with a large range of viscosity. The selected cutoff value of 7 was similar to that of [51].
In the animal study, we observed that the slope of the phase velocity vs. frequency decreased between mid and late pregnancy. Shear wave dispersion in the cervix can be attributed to the layered structure of collagen bundles, and to the presence of the viscous components of the extracellular matrix such as glycosaminoglycans. The collagen bundles are arranged in approximately three bands: an inner band with bundles aligned radially from the canal, a middle one with circumferential bundles, and an outer one with longitudinal bundles parallel to the canal [3]. Each band, with thickness of a few millimeters, can guide the propagation of the shear waves, as observed in tendon [26]. As pregnancy progresses, the layered collagen structure is lost and, therefore, shear wave dispersion caused by guided wave phenomena should decrease.
In addition to the bands of collagen, the thickness of the anterior and posterior portion of the cervix can lead to shear wave dispersion to the presence of Lamb waves. The average halfthickness (thickness of either the anterior or posterior portions) were 2.3 ± 0.3 and 1.1 ± 0.2 cm, respectively. We measured the half thickness instead of the full thickness because each half portion (anterior and posterior) is bounded by the cervical canal and the surrounding fascia. Because of this, we restricted the analysis of shear wave propagation to either the anterior or the posterior portions. To study the correlation between shear wave dispersion and the cervical thickness, we added a plot of the parameter S and cervical thickness ( Figure 11). Even if the  wavelength is similar to the cervix half thickness, we did not find any correlation between the values of the S parameter and the cervical thickness and length (as indicated by the Pearson's correlation coefficient ρ 0.07 (p 0.71)). Also, the concentration of viscous components such as hyaluronan, decorin, and versican varies during pregnancy [52]. Thus, changes in shear wave dispersion during pregnancy result from competing phenomena. Taken together, our results suggest that the structural component of dispersion may have a larger influence on the observed changes between weeks 10 and 23 in the pregnant NHP cervix. To deepen our understanding of the influence of the different structural and molecular components of the cervix on dispersion measurements through SWEI, our laboratory is developing simulation tools based on the finite element methods.
After applying the criteria for reliability of dispersion measurements, eight of the 18 subjects showed reliable estimates in week 10, and 15 in week 23. The larger rejection of data in week 10 can be attributed to the fact that the NHP cervix is stiffer than the human cervix and other organs in which SWEI has been successfully applied, like the liver. The larger stiffness leads to noisier shear wave displacement values. Furthermore, the rejection of a large number of data is common in the initial steps of the development of quantitative imaging. For example, Castañeda-Martinez et al. [53] rejected data from 11 out of 16 subjects while implementing Ultrasound Backscatter Spectroscopy in the neonatal brain. Also, Wang et al. [54] rejected 50% of measurements when using shear wave elasticity imaging in the liver. Because the human cervix is softer than the NHP, we are confident that we will achieve more reliable dispersion estimates in the clinical application of the technique developed in this work.
Various works have reported the use of the shear wave group speed V g as a surrogate for cervical softness in primates and humans [20,23,24,30,[55][56][57][58]. For example [59], reported a decrease in the group speed of shear waves between ripened and unripened ex vivo human cervical samples, 39% for the anterior region and 25% for the posterior region with a p value < 0.05. Also, in a study in patients with a term pregnancy, V g decreased 39% between pre and post-ripening in the context of labor induction, with a p value < 0.001 [18]. Here we observed a reduction in dispersion of 53% with a p value much smaller than the reported changes for the group speed. This indicates that the shear wave dispersion can be more sensitive to cervical remodeling than the group speed. This could be attributed to the more direct relationship between shear wave dispersion and structural properties of the cervix, as well as the optimization of noise reduction strategies.
Using the linear mixed-effects model on the values of S, we found that age and cervix portion (anterior vs. posterior) had a significant influence on shear wave dispersion. Rosado-Mendez et al. [21] observed significant differences in the group speed between the anterior and posterior NHP cervix. Carlson et al. [60] reported similar differences in humans. These differences can be attributed to spatial variations in the microstructural properties of the cervix, a weight gradient created by the fetus, and/or the rotated orientation of the cervix with respect to propagation direction of the shear wave. Another factor could be a pre-compresion of the posterior cervix due to the presence of the transducer. Also, these differences could be attributed to variations in the shear wave frequency content due the attenuation and focusing of the pushing beam. The significance of the mother's age could be attributed to the changes in the production and reception of glycosaminoglycans as the mother age increases [61]. Despite the influence of these factors, the most important effect on shear wave dispersion was the week of gestation, thus indicating the potential of this technique to track cervical remodeling during pregnancy.
This study had various limitations. The noise model that we used in the simulation part was Gaussian and modulated by the maximum shear wave displacement of each displacement field [38]. Although simplistic, this model allowed us to performed a comparative selection of the descriptor of dispersion that maximized the detection of different levels of dispersion. The second limitation is the large fraction of rejected data at week 10. This can be attributed to the stiffer state of the cervix at mid pregnancy compared to late pregnancy, which results in lower signal/noise of displacements estimates. This might be a limitation when performing future comparisons of the early vs. late pregnancy cervix. Since the human cervix is slightly softer than the NHP cervix, we expect that the rejection rate will not be as high in clinical applications of the proposed technique. Finally, the S parameter does not provide a quantitative assessment of the changes in the structural components that determine the dispersive response of the cervix. We are currently exploring the quantification of shear wave attenuation and the fit of rheological models to shear wave dispersion and attenuation to obtain a more detailed characterization of the cervix structural components. Although this strategy has been proposed by other groups, they used simplistic viscoelastic models such as Voigt's [28]. Other models such as the Prony series may provide a more accurate assessment of the structural components.

CONCLUSION
In this study we used shear wave dispersion as a surrogate of cervical remodeling during pregnancy in a non-human primate model. We compared the sensitivity of two dispersion descriptors extracted from either a time-domain or a frequency-domain analysis of shear wave propagation in finite element simulations of viscoelastic media. We observed that the frequency-domain parameter S, defined as the phase velocity slope vs. frequency, in combination with noise-filtering based on singular value decomposition provided better differentiation of viscoelastic materials with different levels of dispersion than other combinations of parameters and filtering strategies. When applied to the NHP model, we observed a significant reduction of shear wave dispersion between mid and late pregnancy in the NHP model, likely due to the disorganization of cervical collagen during pregnancy. Our results suggest the potential of shear wave dispersion as a quantitative imaging biomarker of cervical remodeling. Future studies will apply this analysis to data from longitudinal human studies of normal and preterm pregnancy.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee of the University of Wisconsin-Madison.

AUTHOR CONTRIBUTIONS
AT contributed with data processing, analysis and discussion of results, and manuscript drafting. MP contributed with generation of simulations and discussion of results. HF contributed with data acquisition, analysis and discussion of results, and manuscript drafting. TH contributed to the discussion of results and manuscript drafting. IR-M contributed with data acquisition, supervising data processing, manuscript drafting, and discussion of results.

FUNDING
This work was supported by UNAM-PAPIIT IA104518, IN103219, and IA102320, and by the graduate studies scholarship from the Consejo Nacional de Ciencia y Tecnología to AT. (CVU 863045). This work was also supported by the National Institutes of Health grants T32CA009206 from the National Cancer Institute, and F31HD082911, R21HD061896, R21HD063031, and R01HD072077 from the Eunice Kennedy Shriver National Institute of Child Health and Human Development.