3D ULTRASOUND STRAIN IMAGING OF PUBORECTALIS MUSCLE

—The female pelvic ﬂoor (PF) muscles provide support to the pelvic organs. During delivery, some of these muscles have to stretch up to three times their original length to allow passage of the baby, leading frequently to damage and consequently later-life PF dysfunction (PFD). Three-dimensional (3D) ultrasound (US) imaging can be used to image these muscles and to diagnose the damage by assessing quantitative, geometric and functional information of the muscles through strain imaging. In this study we developed 3D US strain imaging of the PF muscles and explored its application to the puborectalis muscle (PRM), which is one of the major PF muscles. (E-mail: shreya.das@radboudumc.nl) © 2020 The Author(s). Published by Elsevier Inc. on behalf of World Federation for Ultrasound in Medicine & Biology. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


INTRODUCTION
Female pelvic floor (PF) muscles provide support to the pelvic organs by compensating for gravity and abdominal pressure (Hoyte and Damaser, 2016). The set of muscles located in the PF is collectively called the levator ani muscles (LAM). During vaginal delivery, LAM is extended by approximately 245%, allowing the levator hiatus (LH) to widen during crowning (DeLancey et al. 2007;Shek and Dietz 2010;Tubaro et al. 2011;Dixit et al. 2014;Dieter et al. 2015). Vaginal delivery is associated with multiple LAM defects, all of which are risk factors for later-life pelvic floor dysfunction (PFD) (Dalpiaz and Curti 2006;Shek and Dietz 2010;Tubaro et al. 2011;Dixit et al. 2014;Dieter et al. 2015;Notten et al. 2017;de Araujo et al. 2018). PFD comprises disorders that include stress urinary incontinence (SUI), overactive bladder and pelvic organ prolapse (POP) (Bedretdinova et al. 2016;de Araujo et al. 2018). It has been reported that the primary cause of these late-age PFD in women result from the damage of one or more of the LAM, although the symptoms of this damage can manifest years after the actual occurrence (Dietz and Simpson 2008;Dietz 2013;Shek and Dietz 2013).
Imaging plays a crucial role in diagnosis of PFD. Magnetic resonance imaging (MRI) and ultrasound (US) are most frequently used to image the PF. These techniques are mainly used to image the anatomy to diagnose POP or SUI. Segmented organs in MRI or US images are used in biomechanical analysis to gain a better understanding of various pelvic organ disorders or specifically to diagnose POP (Akhondi-Asl et al. 2014;Onal et al. 2014Onal et al. , 2016Nekooeimehr et al. 2016;Wang et al. 2018). Furthermore, biomechanical modeling using imaging as input has been performed to enhance understanding of the anatomy of PF, aid in the diagnosis of various PF disorders and aid in surgical planning for corrective surgeries (Damser et al. 2002;Parikh et al. 2002;Bellemare et al. 2007;Pu et al. 2007;Lee et al. 2009). Diagnosis of POP and SUI is also performed using anatomically significant reference points (e.g., bladder neck and anorectal angle) (Basarab et al. 2011;Czyrnyj et al. 2016). Current studies also include deformation analysis of the PF, which has been performed using PF organs such as the bladder, uterus and rectum (Rahim et al. 2009(Rahim et al. , 2011Ogier et al. 2019). Lastly, research on PF also includes vaginal tactile imaging, magnetic resonance defecography and distribution representation of PF muscles in human motor cortex (Egorov et al. 2010;Costa et al. 2014;Yani et al. 2018). In summary, the PF as a whole has been investigated, yet the PF muscles have not been evaluated individually in detail. As a result, functional information on PF muscles is scarce.
For the detection of PFD, however, it is crucial to have functional information on the PF muscles. Studies that have investigated the PF muscles conclude that the biomechanics of these muscles can be an indicator of PF disorders such as POP (Silva et al. 2017;Hu et al. 2019). The female PF muscles provide support to the pelvic organs. Therefore, apart from understanding the biomechanics of the entire PF, it is also important to determine the function of these muscles when undamaged, to serve as a reference when investigating damaged muscles. Quantification of functional information on and damage to PF muscles might be possible with strain imaging. Strain imaging can be performed using both MRI and US. However, 3D US imaging has several benefits with respect to MRI, for example, its ease of use, portability, minimal discomfort, relatively short period required for acquisition of the data and low price (Shek and Dietz 2010;Tubaro et al. 2011;Shek and Dietz 2013;Dixit et al. 2014;Dieter et al. 2015). Furthermore, good consistency between US and MRI in the imaging of PF muscles has been established (Yan et al. 2017).
US-based strain imaging can be used to investigate muscle movement and has been used extensively to understand and investigate the complex movements of the heart and skeletal muscles (Kalam et al. 2014;Gijsbertse et al. 2017). To date, US-based strain imaging is concentrated predominantly on 2D US images, whereas the PF muscles are complex 3D structures and their movements and deformations constitute an inherent 3D phenomenon. This means that the muscle has to be tracked in three dimensions to accurately quantify the 3D motion and deformation over time. Furthermore, the deformation has to be quantified in multiple directions to fully characterize the deformation of the muscle and to be able to identify dysfunctional or damaged parts of it.
In the work described here, we developed 3D US strain imaging specifically for PF muscles and investigated the puborectalis muscle (PRM), one of the major muscles of the LAM. To the best of our knowledge, no existing study has investigated the function of individual muscles of the female PF. The reason for investigating this muscle is twofold. First, this muscle is frequently damaged during childbirth, and that damage is the primary cause of SUI or POP later in life (Dietz and Simpson 2008;Dietz 2013;Shek and Dietz 2013). Second, the PRM forms the outline of the LH when it is viewed in the top view in US data of the female PF (Grob et al. 2014). The LH is one of the most important parameters assessed in transperineal ultrasound (TPUS) studies. Therefore, studying strain of the PRM first seems to us the most logical choice and would allow us to relate our results to the existing literature.
We hypothesized in this work that it is possible to quantify strain in three dimensions for the deforming PRM using a time series of volumetric US data. Strain imaging was performed in four nulliparous women (n = 4), who were asked to contract (n = 4) their PF muscles or perform both a contraction and a Valsalva maneuver (n = 1). We chose nulliparous women, that is, women who have not yet given birth because these women are considered to have intact PF muscles. Therefore, this study provides insight into how the undamaged muscle strains in three dimensions during activation.
Summarized, the aim of this study was to develop 3D US strain imaging of the PF muscles and to explore its application to the PRM, which is one of the PF muscles.

Data acquisition
Dynamic 3D TPUS volumes were acquired using a Philips X6-1 matrix transducer connected to an EPIQ 7G US machine (Philips Healthcare, Bothell, WA, USA), at the University Medical Centre (UMC), Utrecht, The Netherlands. Total acquisition length was 11À15 s at a volume rate of 1.5 Hz. US data were obtained over time for female volunteers (n = 4) who had never given birth (nulliparous). These women had overactive PFs, which is chronically raised pelvic muscle tone. The PRMs in these women were intact and undamaged, as confirmed by the clinician. US volumes were recorded during two types of exercise. During the first exercise, contraction, the women were asked to actively contract their PF muscles, commencing and ending with the muscles in a state of rest. During the second exercise, they performed a Valsalva maneuver. A Valsalva maneuver is performed by a moderately forceful exhalation against a closed airway. It is used to increase the abdominal pressure, which causes LAM distension and allows the clinician to assess the full extent of a POP (Hoyte and Damaser 2016). Data acquisition commenced with rest and ended at maximum Valsalva maneuver. The data were stored in the Digital Imaging and Communications in Medicine (DICOM) format. Table 1 summarizes the demographic characteristics (age and body mass index) and the exercises performed by the included volunteers.
The Medical Research Ethics Committee of UMC Utrecht exempted the project from approval, and all volunteers signed appropriate research consent forms.

PF imaged through TPUS
Figures 1 and 2 illustrate US data acquired from a PF in the rest state, as imaged with TPUS. In this example, we can observe the PRM in the rest state. The transducer was positioned against the PF while the women were in supine position.
Both the bone pubic symphysis (PS) and the PRM can be visualized with ease in the sagittal view, as can the surrounding PF organs. In Figures 1 and 2, the PRM is the yellow region bordering the anal canal.

Data processing
The block diagram in Figure 3 illustrates the processing steps performed to calculate strain volumes of the PRM. The input comprised the recorded dynamic US volumes and the region of interest (ROI) for which strain was calculated. To obtain the ROI, the PRM was manually segmented by an experienced clinician in the initial US volume (rest before exercise) (van den Noort et al. 2018). The output was a set of accumulated echo volumes as a function of time. The influence of segmentation of the ROI was varied by decreasing the ROI, and displacement estimates were calculated using these volumes. The difference in the values obtained at these different sizes was identical for corresponding voxels. The processing sequence can be divided into four steps: volumetric data preparation, intervolume displacement estimations, tracking (involving an update of the ROI) and strain calculations.
Each processing step is explained in detail hereafter.

Volumetric data preparation
The first of the two inputs for the work was the data from the US machine, which were in DICOM format. These data were first converted to a rectilinear format with ".fld" extension using a proprietary software called QLAB, Version 10.8 (Philips Healthcare, Andover, MA, USA). Conversion of the data was performed to allow import in MATLAB R2018 a (The MathWorks, Inc., Natick, MA, USA), which was the program we used to develop our 3D strain analysis software. The total number of volumes per data set was 22, and each of the 3D volumes contained 352 £ 229 £ 277 (X £ Y £ Z) pixels which were uniformly sampled at distances of 0.42 £ 0.60 £ 0.34 mm (dx £ dy £ dz).

Intervolume displacement estimations
The next step was to calculate intervolumetric displacements. Therefore, displacements were estimated between the first two volumes within the initial ROI (illustrated in Figs. 1 and 2 for volunteer 1). Intervolumetric displacements for each pair of subsequent volumes were estimated with a 3D normalized crosscorrelation algorithm (Gijsbertse et al. 2017;Hendriks et al. 2016) optimized for PF muscles and the US system used in this study.  In this algorithm, two subsequently recorded volumes were subdivided into 3D blocks called kernels and templates. The kernel and template sizes used for those volumes were 111 £ 81 £ 41 and 51 £ 51 £ 11 pixels, respectively. The kernels were matched on the templates, and the locations of the 3D cross-correlation peaks were calculated. These locations of peaks indicated the displacements between the two blocks. To estimate subsample displacements, the cross-correlation peaks were interpolated (parabolic fit) (Hendriks et al. 2016). The displacement estimates were finally filtered using a 3D median filter.

Tracking
As the PRM changes (position and shape) from volume to volume, that is, changes with time, the position and shape of the ROI for displacement estimation also needs to be updated over time. Otherwise, displacement estimation would no longer be performed for PRM tissue only, but would gradually shift to the surrounding tissue. The ROI coordinates (position of the manually segmented PRM) of the next volume were updated using the displacement estimates as calculated in the previous step. In the next step, displacements were calculated between the next two subsequently acquired volumes, and the ROI was updated. The process of estimating intervolumetric displacements and updating ROI is called tracking (Lopata et al., 2009;Lopata et al. 2010). Tracking began after the displacement estimations between the first two subsequently acquired volumes were estimated using the initial ROI.
The input to this processing step was the filtered intervolume displacement estimates from the previous step. Filtering was done using a median filter (2 £ 2 £ 2 cm) to smoothe the displacement estimates and remove outliers (Hansen et al. 2010;Hendriks et al. 2016). This was required for tracking.
Accumulation of the filtered displacement estimates was performed using the equation where accum_disp_estmts = accumulated displacement estimates in z, x or y direction, inter_vol_dis-p_estmts = intervolume displacement estimates in z, x or y direction and n = number of US volume.  The accumulated displacement estimates were the total movement of the muscle up to the (n + 1)th volume or time point. These displacement estimates were obtained in number of US grid points that a certain index had passed. The indexes were updated with these accumulated displacement estimates.
Because the updated indexes are subsample values, displacement estimates were calculated from the eight surrounding sample points, and the displacement estimates of the subsample points were arrived at by linear interpolation (Fig. 4). In this way, the muscle could be tracked throughout its complete deformation cycle. After each update of the ROI, it was checked visually on the respective US volume to ensure whether it was at the same position as the displaced muscle.
Therefore, for this processing step, the inputs were the intervolume displacement estimates, and the outputs were the updated ROIs.

Strain calculations
Strain calculation was the last processing step. Accumulated displacements were calculated by summing intervolumetric displacements up to each time point. The non-filtered intervolumetric displacements were used and filtered with a median filter kernel (1 £ 1 £ 1 cm) before accumulation (Hansen et al. 2010;Hendriks et al. 2016). A kernel size smaller than that in the tracking step was applied to avoid too much smoothing of the displacements, which would result in the absence of a gradient in the strain calculations. As the displacement estimates were filtered using a different kernel, interpolation was again performed, now for the updated indexes. In the next step, the interpolated displacement estimates were accumulated using eqn (1). These accumulated displacement estimates were used to calculate the 3D strain tensor using a 3D least-squares strain estimator (LSQSE) (Kallel and Ophir 1997).
The contraction direction of the PRM is not coaligned with the rectilinear coordinate system. To determine the major or principal component of the strain that is induced in the contraction or Valsalva maneuver of PRM, principal strains were calculated from the individual strain values in the z, x and y directions (Tuttle 2012). As we have observed from the LSQSE strain that strain for contraction is negative and strain for Valsalva maneuver is positive, we chose the largest negative principal strain component for data during contraction and the largest positive strain component for data during the Valsalva maneuver.

RESULTS
Accumulated displacement estimates are illustrated in Figures 5 and 6, and principal strain results, in Figures  7 and 8, for two of four volunteers. These time points are, respectively, muscle at rest, muscle at maximum contraction and muscle at rest post-contraction, for volunteer 1. In the case of volunteer 4, the time points are rest and maximum Valsalva maneuver. The principal strain magnitudes and principal strain directions are illustrated in the figures.

Accumulated displacement estimates
In volunteer 1, at the time point at which the muscle is at rest (Figure 5a, 5d, 5g, first column), the estimated displacements between the first two volumes are quite low in all directions, which is expected at rest. We observed this for displacement estimates of all volunteers.
In Figure 5b, 5e, 5h (second column) are the accumulated displacement results for maximum contraction. The displacement estimates in the z-direction are the highest, followed by displacement estimates in the ydirection. In the x-direction, we see that displacement estimates are almost zero.
During contraction, in the z-direction, negative displacement estimates mean that the PRM is moving toward the bone PS and, thus, toward the US transducer. In the y-direction, displacement estimates are positive, which means that in this direction, the muscle is moving away from the US transducer. There is very little lateral or side-to-side movement of the muscle, that is, in the xdirection. These movements or lack of movement with respect to the muscle at rest are illustrated in Figure 9bÀd. In Figure 5c, 5f, 5g (third column), we see that the muscle is almost back at rest, and so the accumulated displacement estimates are again back to approximately zero values. The muscle does not completely return to the rest position, because these data sets were acquired in women who have overactive PFs. Thus, these women might take longer to return to the rest position post-contraction.
In the data sets acquired during the Valsalva maneuver, data acquisition was stopped when the muscle reached maximum Valsalva maneuver. The accumulated displacement results are illustrated in Figure 6. We observe that displacements are initially close to zero, as the muscle is at rest.
During the maximum Valsalva maneuver the accumulated displacements are predominantly positive in the z-direction and negative in the x-and y-directions. This indicates elongation of the muscle in Valsalva maneuver as opposed to shortening during contraction. This deformation with respect to the muscle at rest is illustrated in Figure 10bÀd.
The movements of the PRM during contraction and Valsalva maneuver, axial, sagittal and coronal views, are illustrated in the supplementary videos in the Supplementary Data (online only). In these videos, the dark gray area represents the position of the muscle during rest, and the yellow area, the muscle during contraction/ Valsalva maneuver. We can observe that the movement during contraction with respect with the bone PS is complementary in direction compared with that of Valsalva maneuver.

Principal strain values
During contraction, as illustrated in Figure 7, it is observed that the major principal strain becomes more negative with increasing contraction. As the muscle returns to the rest position, the negative strain decreases but does not become zero.
The principal strain component directions change for all volunteers when the muscle contracts from rest and becomes predominantly aligned with muscle fiber direction at maximum contraction. In Figure 7e, 7f, it can be seen that the direction of strain further changes when the muscle returns to the rest state after contraction.
As illustrated in Figure 8, the data for volunteer 4 contain a Valsalva maneuver. In this case, rest remains the same as the rest during contraction, whereas during maximum Valsalva maneuver, strain is positive with a peak value of 60% strain.
The bottom row of Figure 8 illustrates the principal strain component directions. Once again it is observed that the directions change when the muscle changes from rest to maximum Valsalva maneuver. Table 2 lists the spatial means of the principal strain (%) values over the PRM for all data sets. Mean principal strain (%) values in the rest position for all five volunteers were less than 3%. The principal strains at maximum contraction ranged between À8.9% and À41.5%. For the data set for the Valsalva maneuver, the maximum principal strain was 38.6%. At the last time point, namely, rest post-contraction, strain had decreased with respect to earlier time points before but was not equal to that before contraction.

DISCUSSION AND SUMMARY
To our knowledge, this is the first study in which 3D displacement and strain were estimated in the PRM. Normally, when clinicians examine the PF with TPUS, they can only visually examine the motion of the PF and measure the (relative) motion for certain specific anatomic landmarks. In other words, a qualitative assessment can be made. As can be observed, the proposed algorithm allows quantitative determination of strain, locally within the PRM.
We observe from the obtained results that the strain during contraction and the strain during Valsalva maneuver are complementary, which allows the algorithm to distinguish between these two opposite movements of the muscle. We also observe that there are changes in the deformation of the muscle. These deformations are different for contraction compared with Valsalva maneuver. Because of the short data acquisition time, we can observe the muscle returning to almost rest but not complete rest post-contraction. For this initial study, we focused on strain estimation in the undamaged and intact PRM of nulliparous women undergoing voluntary contraction or Valsalva maneuver before focusing on strain estimation in patients with a complex pathology, for example, avulsions.
We observe in the last column of displacement estimates in Figure 5 that the muscle does not return to the exact rest position post-contraction. There are three possible explanations for this: a clinical one, an explanation related to the hardware and a technical explanation. The clinical explanation is that the volunteers from whom the data were acquired had overactive PFs. Therefore, it might be that the PRM will take more time to return to its rest position post-contraction. For example, in volunteer 1, we find the maximum contraction is at volumes 13 and 14, which means that 7 of the 11 s of total data acquisition time had already passed. The PRM does not return to rest within the remaining approximately 4 s. As the volunteers were supine and contracting their PF muscles only during data acquisition, motion of the volunteer or global motion can be ignored. This can also mean a hardware limitation; the time for data acquisition was too short and ended before the muscle had returned to its rest position. Lastly, the technical explanation is that tracking might not be ideal. Small inaccuracies in the displacement estimates accumulate over time to introduce error in tracking.
The complementary sign is visible in the displacement estimates and strain during Valsalva maneuver (Figs. 7 and 8), compared with those observed during contraction (Figs. 5 and 6). It can be observed that the muscle has moved away from the US transducer in the zdirection, and there is clearly a change in shape of the muscle. In the x-direction, we can observe that one end of the muscle has moved more than the other end. This differs from the results for contraction, where there is little or no movement in the x-direction. This might be because the contracting muscle is expected to move toward the bone PS (in effect the US transducer) to which it is attached. It is not expected to move laterally or from side to side in contraction. In the Valsalva maneuver, while the muscle is elongating, we observe that deformation is occurring in all three z, x and y directions. In this volunteer, in the x-direction, one arm of the muscle is manifesting more displacement than the other. In the y-direction, we observe that the muscle has moved toward the transducer. Because we studied the Valsalva maneuver in only one volunteer, these observations cannot be generalized. The presence of these trends needs to be investigated in a large sample size. The PRM is almost uniformly strained (Fig. 6aÀc, first row), indicating contraction, whereas it manifests nonuniform strain when it returns to rest after contraction. A possible explanation might be that in the case of an overactive PF, different parts of the muscle take longer to return to zero strain. As the sample size was small in this study, it should be extended in future studies to investigate whether this trend is present in a large group of women. Also, there is a large variation between the strain (%) values obtained during contraction from the four volunteers in contraction, as illustrated in Table 2. A possible reason could be that different women have different levels of control over their PRMs during contraction. Additionally, deformation of the muscle during contraction might also vary per woman. Larger sample sizes in future studies could provide a range of strain from minimum to maximum for undamaged PRM.
We can observe the directions of the principal strains in Figures 7 and 8dÀf, second row. It can be observed that the strains are in the direction of muscle fiber orientation.

Clinical significance
First, knowledge of the exact length by which the muscle has moved and deformed is beneficial in that we can assess numerically how far the woman can move her muscle voluntarily. Thereafter, we can compare how much movement is expected in a normal undamaged muscle compared with a damaged muscle. Moreover, it would provide clinicians and pelvic physiotherapists with a quantitative tool to follow patients' improvements during treatment (e.g., PF muscle training). Second, when the muscle is observed in the image displayed in the US machine, it is difficult to assess quantitatively which part of the muscle is displaced more and which is displaced less. As illustrated in the results, for undamaged muscles, it is possible to know which part of the muscle is displaced more from the different colors of the displacement estimates in the figures. For example, in the case of displacement estimates in the z-direction, for volunteer 1 (Fig. 5e), the "sling" of the PRM moved more than the two ends that are attached to the bone PS. In the results for volunteer 4 (Fig. 6), in the Valsalva maneuver, we observe that different parts of the PRM move dissimilarly in all three z, x and y directions. These observations of the figures can give us an idea about how the undamaged muscle moves.   Thereafter, a comparison can be made with damaged muscle movement. Lastly is the exact time point, or more specifically the volume, at which the muscle, for certain data at maximum contraction/Valsalva maneuver, can be determined. It can be useful in TPUS as a means of automatically arriving at the specific volume for maximum contraction/Valsalva maneuver, thus reducing a source of variability in TPUS assessments.
In this study, the primary reason for calculating strain induced in the PRM was to quantify the deformation and strain in undamaged PRM. When there is damage or scar tissue formation in the muscle, it might be of clinical significance to assess the exact position of the damage through the different strain (%) values in different parts of the muscle along with the directions of the strain values. Thereafter, it might also be possible to quantify in three dimensions which part of the muscle is damaged.
To further investigate PF muscles through 3D strain, future studies should include larger sample sizes for both undamaged PRM and complex PRM pathologies. Other LAMs should also be investigated to learn how the muscles behave in relation to each other.