The cervix is essential to a healthy pregnancy as it must bear the increasing load caused by the growing fetus. Preterm birth is suspected to be caused by the premature softening and mechanical failure of the cervix. The objective of this paper is to measure the anisotropic mechanical properties of human cervical tissue using indentation and video extensometry. The human cervix is a layered structure, where its thick stromal core contains preferentially aligned collagen fibers embedded in a soft ground substance. The fiber composite nature of the tissue provides resistance to the complex three-dimensional loading environment of pregnancy. In this work, we detail an indentation mechanical test to obtain the force and deformation response during loading which closely matches in vivo conditions. We postulate a constitutive material model to describe the equilibrium material behavior to ramp-hold indentation, and we use an inverse finite element method based on genetic algorithm (GA) optimization to determine best-fit material parameters. We report the material properties of human cervical slices taken at different anatomical locations from women of different obstetric backgrounds. In this cohort of patients, the anterior internal os (the area where the cervix meets the uterus) of the cervix is stiffer than the anterior external os (the area closest to the vagina). The anatomic anterior and posterior quadrants of cervical tissue are more anisotropic than the left and right quadrants. There is no significant difference in material properties between samples of different parities (number of pregnancies reaching viable gestation age).

Introduction

The cervix, the soft tissue structure at the bottom of the uterus, plays a mechanical role in achieving a normal, term pregnancy [1,2]. For the majority of gestation, the cervix must remain closed and strong enough to hold the increasing mechanical load of the growing fetus. At the time of parturition, it must be compliant enough to deform and dilate to allow for delivery of the fetus. To achieve compliance, the cervix undergoes a material remodeling process not completely understood, especially for human tissue.

Deficient material properties and premature cervical remodeling are hypothesized to result in spontaneous preterm birth [3], a leading cause of neonatal morbidity and neonatal death [4]. Preterm birth is a major global health issue with an estimated 15 × 106 babies born preterm around the world [4], where a majority of these are spontaneous [3]. It is vital we obtain a better understanding of the material characteristics of the human cervix to discern its complex mechanical function during pregnancy and to identify mechanisms of normal and abnormal remodeling. In this study, we present experimental results to characterize the passive, anisotropic material behavior of the human cervix to indentation loading. The aim here is to establish a baseline understanding of how the collagen fiber architecture contributes to the compressive material response, using nonpregnant (NP) cervical specimens of various parities. Future studies will seek to characterize pregnant tissue samples from patients with and without a history of preterm birth with hopes to understand the mechanical mechanisms of premature softening.

The human cervix can be idealized as a thick-walled cylinder. The overall size of the cervix has been reported to be around 3 cm long and 2.5 cm in diameter [5]. These dimensions, along with the size of the inner canal, vary widely between patients [6]. The overall range of cervical volume and size for a large cohort of women with varied demographics and obstetric backgrounds is not known.

The cervix is a layered tissue, where the outer edges are lined with epithelial cells and the inner canal is lined with a thin layer of mucosal cells. The dense core of the cervix is called the stroma and it provides the overall mechanical integrity. The tissue in the cervical stroma is a dense network of collagen fibers (including but likely not limited to types I and III collagen fibers) composing 34–77% of the dry weight [7,8]. This large range in collagen values may have something to do with regional differences in the tissue recently found through immunohistochemistry staining of smooth muscle cells [9]. The internal os contains a significant amount of smooth muscle (50–60%) dispersed throughout the stroma that is circumferentially oriented around the endocervical canal. The content of smooth muscle gradually decreases to 10% as one moves toward the external os [9].

The collagen fibers form a hierarchical network embedded in a viscous elastic ground substance of negatively charged glycosaminoglycans and other proteins [1]. Similar to smooth muscle in the cervix, the collagen fiber network is reported to be anisotropic with different preferred orientations in distinct anatomical regions within the cervix. According to previous studies using optical coherence tomography (OCT) [6], X-ray diffraction [10], and magnetic resonance imaging [11], the inner and outermost layers of the cervix contain fibers aligned parallel to the inner canal and the middle layer contains fibers aligned circumferentially around the inner canal. These regions are overlapping, with seamless transitions, creating a complex architecture of dispersed fibers that are able to withstand load in multiple directions [12]. As seen in the most recent of these studies [6], circumferential fibers occupy both the middle and outer layers of the cervix, covering the majority of the stroma. In the mode of deformation induced by spherical indentation directed parallel to the inner canal, as done in this study, the circumferential fibers engage to provide compressive support [13].

It is the hypothesized cervical softening occurs by altering the extracellular matrix (ECM) content and its multiscale organization [1]. Our group has previously published uni-axial tension and compression [13,14] and indentation [15] mechanical testing data on the human cervix. These studies report the mechanical properties of cervical tissue are anisotropic, time-dependent, tension/compression asymmetric, and influenced by the number of times a woman was pregnant. Moreover, pregnant tissue is considerably softer and more compliant than nonpregnant tissue. In equilibrium tension, at strain levels up to 30%, the material parameter associated with the fiber stiffness decreases by two orders of magnitude for pregnant tissue (∼100 kPa to ∼ 1 kPa) [13,16]. In equilibrium compression, the differences between nonpregnant and pregnant tissue are more subtle given its low compressive stiffness (∼1 kPa) [13,15,16]. These initial mechanical tests provide a starting framework to understand the material behavior of human cervical tissue. A more advanced understanding of the collagen fiber architecture [6,10,11] and the complex loading conditions in pregnancy [12,17] motivates the development of a mechanical test to characterize tissue anisotropy in large compressive deformation.

Indentation is a popular method to mechanically test soft tissue [1823]. Compared with tension or other mechanical test methodologies, samples need minimal preparation and often indentation is nondestructive. Additionally, with the minimal sample preparation, samples are able to preserve features of the collagen fiber architecture. On the other hand, due to the complex boundary conditions of indentation, it is challenging to extract a material stress–strain response from the experimental data. Hence, in order to quantify material parameters from indentation, this study presents an experimental and computational method to generate material properties for an anisotropic fibrous tissue.

Specifically, the objective of this paper is to quantify the material properties of nonpregnant human cervical tissue samples using information about its collagen fiber architecture [6] and spherical indentation. The focus of this paper is to characterize the anisotropic material behavior at various anatomical locations of the cervix after the transient response of the tissue has subsided. Here, we present an indentation test coupled with digital image correlation (DIC) to capture two-dimensional (2D) strain responses. Using the indentation force–displacement, 2D strain field, and collagen fiber architecture OCT data [6], inverse finite element analysis (IFEA) based on genetic algorithm (GA) optimization is used to quantify material parameters of a fiber composite material. GA is widely proven to have a better global optimization ability compared with other algorithms [2427]. Our results show human cervical tissue exhibits anisotropic and time-dependent material behavior. Further, cervical tissue at the anatomic anterior of the internal os is significantly stiffer than cervical tissue at the anatomic anterior of the external os. The anterior and posterior quadrants have tighter aligned collagen fibers than the left and right quadrants. In this small cohort of patients, there is no significant difference found between cervical samples from women with different parities.

Method

Sample Preparation.

Cervical tissue samples were collected from seven healthy NP women undergoing hysterectomy for benign indications using an Institutional Review Boards approved protocol at Columbia University Irving Medical Center. The demographic information for each patient is listed in Table 1. Immediately after hysterectomy, the cervix was separated from the uterus at the area of the internal os (defined as where the uterine arteries meet the cervix/uterus) using a custom slicer [15] (Fig. 1).

Transverse cervical tissue slices were then obtained from the area of the internal and external os using an established protocol [6]. Tissue slices were immediately frozen on dry ice and stored at −80 °C until they were used for mechanical testing. Due to the difficulty of slicing, the thicknesses of the slices varied from 1.76–7.07 mm (the average was 3.82 mm and the standard variance was 1.29 mm). Slices were numbered, starting at the internal os (i.e., the top of the cervix closest to the uterus). Depending on the dimension and integrity of different cervical slices, 5–10 slices were collected. Two to three slices were used for the indentation tests and the remaining were used for other experiments. A more detailed protocol used for sample collection and preparation was described in our earlier work [6].

Mechanical Tests.

Twenty-four hours before imaging, each slice was defrosted in 4 °C and allowed to equilibrate in phosphate buffered saline (PBS). The surface was microtomed, fiber architecture measured via OCT [6], and then restored at −80 °C. At the time of the mechanical test, the slice was equilibrated in PBS at 4 °C for 12 h. Each slice was then weighed and measured. The bottom surface of the cervical slice specimen was dried by KimWipes and speckled with Verhoeff's elastic stain (VEG) by an airbrush (Harder & Steenbeck) with a 0.2 mm tip (Figs. 2(a)2(c)).

After speckling, each specimen was placed in a petri dish with glue applied to the outer edge of the specimen to ensure the specimen did not float. The petri dish was then filled with PBS and the specimen was equilibrated for at least half an hour at room temperature. The PBS bath was placed on the rigid surface of a universal testing machine (Instron, Inc. Microtester 5948, Norwood, MA) equipped with a 6 mm diameter sphere indenter attached to a 5 N load cell (Instron, Inc., Norwood, MA, 0.25% accuracy of indicated load).

The zero point for indentation contact was found following our previous protocol [15]. Once the slice was recovered, a displacement-controlled, multilevel ramp-hold indentation test was conducted in PBS. There were a total of four ramp-hold levels with indentation depths determined by percentages of the specimen thickness at the indentation location (Fig. 2(b)). For thickness more than 5 mm, the indentation depths were 12.5%, 25%, 37.5%, and 50% of the thickness. For thickness less than 5 mm, the indentation depths were 15%, 30%, 45%, and 60% of the thickness, in order to achieve a more significant change in force and strain response between levels. The ramping rate was always 1% thickness per second. The hold times were 480 s, 600 s, 720 s, 960 s for the four levels, respectively (Fig. 3(a)). Load-time (N·s) data (Fig. 3(c)) were recorded using a material tester software (Instron, Inc., Blue Hill version 3.11.1209). Each indentation test was performed in the midstroma area in the center of each quadrant (posterior (P), anterior (A), left (L), and right (R)) of each slice, sufficiently away from the glued boundaries. The midstromal region was defined as the region apart from the fascia at the outer area and the mucosa around the inner canal (Fig. 2(b)).

Digital Image Correlation.

To measure the deformation of the bottom surface of the slice (Fig. 2(a)), the slice was imaged during the indentation test by a charge-coupled device camera (Point Gray Grasshopper, GRAS-50S5M-C75 mm, f/4 lens) using Vic-Snap (Correlation Solutions, v2010, Irmo, SC). The lowest possible aperture, producing high contrast images, was used to maximize the depth of view. The PBS bath was placed on a rig with a transparent window. The light was cast by a fiber optic illuminator (Cole-Parmer, model 41500-50) from the top of the rig. A 90-deg prism was fixed below the window to image the bottom surface of the specimen from the front. During ramps, the acquisition rate was two images per second. During holds, the acquisition rate was one image per 10 s.

The images were processed for digital image correlation (DIC) by Vic-2D (Correlated Solutions, v6) with the incremental method. In the incremental correlation mode, each image is correlated with the previous one. The area of interest was selected around the center of indentation. A subset size was selected for each specimen, between 29 and 59 pixels, depending on the speckle pattern such that within each subset there were enough unique features for correlation. Step size, which controls the spacing of the points analyzed during correlation, was set to one-quarter of the subset size as suggested by the software manual. A starting point was selected in an area without much deformation so the speckle pattern was preserved throughout the experiments. The displacement field in mm was converted from the pixel field. The conversion factor was calculated using the calibration picture of a standard ruler, which was slightly different for each case. The average conversion was 0.017, i.e., δmm=0.017δpix. After DIC, Lagrange strains in horizontal (exxDIC) and vertical (eyyDIC) directions and shear strain (exyDIC) were calculated using 90% centered Gaussian filter and extracted (Fig. 4(a)). The x- and y-directions indicated here are in the experimental frame, and all of the experimental DIC 2D strain fields (Fig. 4(a)) are reported in this frame. The relationship between the experimental frame and the frame used for IFEA is given in the Appendix.

Finite Element Analysis Model.

For the finite element analysis (FEA) model, each specimen was represented as a half-disk, and the indenter was represented as one-quarter of a 6 mm—diameter hollow sphere (Fig. 4(b)). Symmetric boundary conditions and the indentation profiles were set up according to the real experimental conditions with appropriate simplification. FEA was run using febio.2

Model of Cervical Slice.

Sample geometry was taken directly from the charge-coupled device camera images. Figure 4(a) describes the measurements of the following critical points:

  • (a)

    one end point of the canal;

  • (b)

    the other endpoint of the canal;

  • (c)

    midpoint of AB;

  • (d)

    indenter spot;

  • (e)

    glued boundary edge;

  • (f)

    sample edge.

Sample dimensions were converted into a solid model representation using a custom matlab (v2016b) script. The representative geometry of the sample was created as a flat cylinder hollow in the middle to represent the cervical slice with a canal. Figure 4(b) describes the solid model using the following parameters:

  • (a) effective radius of the inner canal;

  • (b1) distance between the indenter spot and the edge of the canal;

  • (b2) distance between the indenter spot and the glued boundary;

  • (c) distance between the glued boundary and the edge of the sample.

The solid model of the cervical slice was meshed using HEX-8 elements. HEX-8 elements have comparable accuracy with TET-10 quadrant elements but are computationally cheaper [28]. The mesh was divided into fine and coarse areas. The fine area was defined as the area with the center as the indentation point D and the refinement radii Rfc,Rfr, and Rfl in the circumferential, radial, and longitudinal direction, respectively. In this work, we used Rfc=3mm,Rfr=3mm, and Rfl=1.2mm, respectively, which were accurate enough based on mesh refinement analysis (not shown here). The other area was defined as the coarse area, and the nodes were equally distributed for both the coarse and fine areas. For the coarse area, the distance between two neighboring nodes was < 1 mm while for the fine area the distance between two neighboring nodes was < 0.4 mm. The number of elements ranged from 1680 to 11,583 depending on specimen dimensions. Using these criteria and geometric parameters, the elements of the fine area were quasi-cubes with the side-length = 0.4 mm. The mesh parameters produced precise enough simulations, with acceptable computation times.

Model of Indenter.

The indenter was meshed using a butterfly distribution of a half hollow sphere with thickness = 0.1 mm, where the size of the elements in the contact area was about 0.4 mm × 0.4 mm. The size of the elements at the center of the top surface of the specimen was about 17% of the radius of the indenter. The indenter was set as a rigid body.

Contact, Boundary Conditions, and Loading Profile.

The contact was selected as “sliding–tension–compression” between the surface of the indenter and the surface of the sample, which represented a frictionless interface. For the bottom surface, z displacement of the ranges b1 and b2 was constrained and all the displacements of the range c were constrained. The sample and indenter geometries, the boundary conditions, and the loading pattern were symmetric to the plane indicated in Fig. 4(b). Hence, a zero-displacement boundary condition normal to the plane was applied to save on computation time.

The displacement loading profile was set to the exact experimental conditions as stated above. The process was divided into four equal steps to represent the four ramp-hold levels, where hold times were adequate for the equilibrium analysis.

Constitutive Model.

The cervical stroma is modeled as an anisotropic fiber composite material, where a continuously distributed fiber network is embedded in a compressible neo-Hookean ground substance. A version of this model for cervical tissue was previously published [29]. Briefly, the total Helmholtz free energy density Ψ is the additive contribution from the fibrous collagen network ΨCOL and a compressible neo-Hookean ground substance ΨGS of noncollagenous ECM components, which is expressed by 
Ψ(F,n,b)=ΨGS(F)+ΨCOL(F,n,b)
(1)
where F is the deformation tensor, and n is the preferential direction of the fiber distribution and b is the von-Mises concentration parameter (b >0) of the fiber.

Each sample in this study was imaged using OCT before mechanical tests to determine the collagen fiber orientation n in the plane perpendicular to the inner canal [6], and the dispersion parameter b was left as a fitting parameter. The OCT maps reveal the stroma region of the human cervix has preferentially aligned fibers in the circumferential direction. These circumferential fibers are a dominant feature, occupying the majority of the stroma region. A few specimens also have a small region of radial fibers immediately adjacent to the inner canal. These radial fibers are not modeled here because of already existing complexities at the boundary. Additionally, we are not considering axial fibers near the boundaries, as found in other studies [10,11], because of similar boundary considerations and axial fibers do not participate in holding the compressive load for this particular indentation deformation. Hence, our analysis will exclusively model a single family of circumferential fibers dispersed about its preferential direction and the material parameter optimization will focus on regions away from boundary edges.

The strain energy function for the neo-Hookean ground substance ΨGS is 
ΨGS=μ2(I13)μlnJ+λ2(lnJ2)
(2)
where μ and λ are Lamé parameters, which could be converted to the more frequently used Young's modulus E and Poisson's ratio ν. I1 is the first invariant of right Cauchy–Green stretch tensor C=FTF, and J=detF is the relative volume.
The fiber network is modeled as a preferentially aligned single family of fibers circling the inner canal [30], with a three-dimensional von-Mises fiber distribution [6]. In this model, fibers can only support tensile loading. The free energy density of the collagen network is given by 
ΨCOL(F,n,b)=02π0πH(In1)ΨnCOL(In)eb(2n21)2πI0(b)sinϕdϕdθ
(3)
where H(⋅) is the Heaviside unit step function, and In=n·C·n. ϕ and θ are the spherical coordinates. The detailed definition of the modified zeroth-order Bessel function I0 is found in our previous work [6].
The strain density function of a single fiber is given by the exponential form 
ΨnCOL=ξ2α[eα(In1)21]
(4)
where the stress-like parameter ξ > 0 and the exponential index α0.

Based on this constitutive model, five material parameters E, ν, b, ξ, and α are used to characterize the material behavior. E is related to the stiffness of the ground substance, ν is related to the Poisson's effect of the ground substance, ξ represents a fiber stiffness, b represents fiber dispersion, and α controls the shape of the exponential.

Inverse Finite Element Analysis With Genetic Algorithm.

The material model was fit to both the force–displacement data and the equilibrium normal and shear strains for an average of points within 0.1 mm of the point directly under the indenter tip (The error bar for the strain in Fig. 5 represents the spread of these data). The five material parameters (E, ν, ξ, b, and α) were optimized, with the following boundaries listed Table 2. Material parameter boundary values were determined from previous uni-axial tension and compression experiments of human cervical tissue [16]. This fit was done for each quadrant indention location around the cervical sample (e.g., anterior, posterior, left, and right as indicated in Fig. 2(b)). Hence, the experimental DIC strain measurements were rotated to match the coordinate frame of the FEA when needed. The detailed method to transform DIC strain measurements to the FEA coordinate frame is detailed in the Appendix. In this way, the transformed strain measurements (exxEXP and eyyEXP) coincided with the Lagrangian strains in circumferential and radial direction of the tissue, respectively, at the spot right under the indenter. Best-fit material parameters [E, ν,ξ, b, and α] were found by optimizing the following objective function 
Ξ(E,ν,ξ,b,α)=i=1N|exxFEAiexxEXPienormEXPN|+i=1N|eyyFEAieyyEXPienormEXPN|+i=1N|FFEAiFEXPiFEXPN|
(5)

Here, the superscript i represents the number of levels, N was the total number of the levels (which is four in this case), and F represented the equilibrium force on the indenter. enormEXPN=[exxEXPN]2+[eyyEXPN]2 was the normalized strain factor of the last level, used to normalize the strain error. The terms with the superscript “FEA” and “EXP” denoted the strains and force response of the FEA and the experiment, respectively.

A GA using a custom matlab script was implemented to optimize Ξ. All computations were run on Columbia University's high-performance computers (Habanero).3 The whole process of each iteration of GA was divided into three operations: selection, crossover, and mutation. The GA started from the zeroth generation G0: a randomly generated population P0(p10,p20,p30,,pn0) in the search space. In this problem, the size of the population was chosen as n =4, which proved to have enough stability and accuracy. An individualpj0 in this population meant a randomly generated combination of parameters (E, ν, ξ, b, and α), which was therefore a vector ((pj0)k,k=1,,5). E and ξ were chosen in the logarithmic scale to avoid them from falling into the dilemma of finding only from the greatest order because the ranges of those parameters were always beyond several orders. In each generation, there were two steps, crossover and mutation. In the crossover step, each parameter from a parent had 50% chance to be inherited to the next generation. In the mutation step, the parents' parameters were modified by a random value, which obeys a normal distribution (defined as the mutation rate). Six children were generated from four parents (C42). Ten individuals, including all parents and children, then went through a tournament selection process in which the individual with the lowest fitness score (Ξ in Eq. (5)) has higher chance to survive and exactly four individuals would survive and become parents for the next generation. A total of 40 generations was iterated. The initial mutation factor μm = 0.15. The mutation factor was divided by 1.3 when the best fitness did not improve for more than five generations. The GA process was validated by comparing the results with the exhaustive method for a randomly chosen case and showed a stable and accurate result (error less than 5%).

Model Validation and Model Error.

In this study, force and strain data from directly under the indenter were used for the IFEA optimization. In order to validate the predictive capability of the model, experimental strain measurements not used for the IFEA fit were compared to the calculated FEA results using the optimized material parameters. In other words, all of the strain measurements except for the point under the indenter tip was used for model validation. For each element of the bottom surface of the FEA model, the calculated strains (exxFEA,eyyFEA, and exyFEA) were compared with the strains of the corresponding position of the DIC results (exxEXP,eyyEXP, and exyEXP). The error for each element of FEA was calculated using Err=|(eFEAeEXP)/eEXP|.

Statistical Analysis.

A group of analysis of variance (one-way analysis of variance) was performed in matlab and the optimized material parameters (E,ν,ξ,b,α) were compared among all the samples with several different criteria. The samples were divided into different groups according to (1) corresponding patient's parity number (parity < 2: Ncervix = 4, parity > 2: Ncervix = 3, where Ncervix is the number of patients involved in the statistics), (2) the anatomical quadrants (anterior and posterior N = 26, left and right N = 26, where N represents the number of indentation locations if not explicitly expressed.), and the location at the cervix (external os Nslice = 7, internal os Nslice = 6, and Nslice is the number of slices used in statistics.), and (3) the samples of the internal os and the external os were also compared within specific quadrants.

Among all the samples available for the indentation tests, seven slices were from four patients who delivered fewer than 2 times, and six slices were from three patients who delivered more than two times. Patients delivered 0 or 1 time were defined as patients with fewer pregnancies, and patients delivered more than 2 times were defined as patients with more pregnancies. Coincidentally, we did not have samples from patients who delivered exactly 2 times. Therefore, the two groups were specifically divided with regard to the different number of parities. The slices with serial number less than half of the total number of the slices were viewed as slices from the internal os, otherwise, the slices were viewed as slices from the external os.

Results

Indentation Force–Strain Response.

Cervical tissue displays a force–relaxation response to ramp-holds in indentation displacement (Fig. 6), and the indentation force response has a nonlinear force–displacement equilibrium response (Figs. 3 and 5(b)).

Sample geometry of the axial cervical slices are circular (Fig. 2(b)), with indentation spots located such that the tissue horizontal or vertical strain measurements under the indenter correspond either to the circumferential or radial directions. However, often times axial slices are slightly misshapen and the DIC coordinate frame did not line up with an anatomical direction (a limitation of this study discussed in Sec. 4.5). Therefore, examples of the measured 2D strain field are given in the DIC-frame (Fig. 4(a)). The 2D strain fields (exxDIC,eyyDIC, and exyDIC) have a complex pattern governed by the unique shape of the human cervix, the boundary conditions of the experimental setup, and the preferred direction of the collagen fibers (Fig. 6). Finally, regions of the 2D strain map indicate compressive strains in both the exxDIC and eyyDIC. These zones of compression are most likely due to the glued boundary conditions.

Focusing on the midstromal region, away from the boundaries, exxDIC and eyyDIC are at a maximum under the indenter (Fig. 6). Whereas, in this region, the exyDIC is at a minimum. After rotating the DIC strain measurements to the FEA-frame such that the indenter location is aligned with an anatomical direction (see the Appendix), the magnitude of the maximum eyyEXP is much higher than the maximum of exxEXP (Fig. 5(a)), indicating a degree of anisotropy for the cervix. For the region under the indenter, exxEXP and eyyEXP represent the strain response in the circumferential and radial directions, respectively. Hence, the higher eyyEXP gives evidence of an anisotropic material behavior where the main fiber direction is in the circumferential direction, matching our previously published OCT results [6].

Inverse Finite Element Analysis Results: Best-Fit Material Properties.

The optimized and averaged material parameters for the human cervical tissue specimens are listed in Table 3. Additionally, all optimized material parameters for every indentation location for every cervical slice is in Supplementary Data on the ASME Digital Collection. Values in the table reflect the average of all data (52 indentation locations in total, where 4 quadrant locations per slice, 13 different slices from the internal and external os, and from 7 patients as listed in Table 1). All corresponding raw equilibrium force, displacement, and DIC strain data used for the fitting optimization are archived online with Columbia Academic Commons.4

The GA optimization for all the samples converged. For each case, three rounds of GA were conducted and the results show the GA produces stable solutions for this problem. The stiffness parameter E and ξ for all samples are within the range of [10−1, 101] kPa, except for five cases with 10 < ξ < 33.2 kPa. The average Poisson's ratio ν is 0.22, indicating the compressibility of the cervical samples. The fiber concentration parameter b varies sample to sample and quadrant to quadrant. The average parameter α is 5.7 indicating a level of fiber nonlinearity.

The anisotropic material model did a good job at describing the experimental 2D normal strain field and indentation force data after the transient responses have died away (Fig. 5). The errors calculated by the objective function average Ξ = 0.58 ± 0.22, and 51 out of 52 cases are with Ξ < 1. As a validation, the features of the normal 2D strain field (exx and eyy) of the bottom surface are captured well by the calculated FEA strain fields (Fig. 7). However, the shear strain field is not captured well by the computational model (limitations are discussed in Sec. 4.5).

Comparing Samples.

The anterior quadrant of the cervical samples of the internal os is significantly stiffer than the anterior quadrant of the cervical samples of the external os. Both the parameter E and α have a significant difference between the internal os and the external os, with the p values equal to 0.0037 and 0.0019, respectively (Figs. 8(a) and 8(b)). No other mechanical parameters are found a significant difference between the internal and the external os.

The anterior and posterior quadrants have a significantly (p = 0.0185) higher concentration parameter b compared to the left and right quadrants (Fig. 8(c)). A lower concentration parameter is associated with a higher fiber dispersion. In other words, the anterior and posterior quadrants have more preferentially aligned fibers. No other mechanical parameters are found to have a significant difference between the anterior and posterior quadrants and the left and right quadrants.

There is no significant difference found in any material parameters between the groups of cervical samples from women with different parity numbers. However, considering the small patient numbers, this study is not statistically powered enough to make a definite conclusion on how parity affects material parameters. This is discussed in Limitations (Sec. 4.5).

Discussion

In this study, we develop a spherical indentation mechanical test to measure the anisotropic material behavior of nonpregnant human cervical tissue. We improve on previous indentation methods by: (1) including measurements of the 2D strain field of the sample against a rigid substrate to reveal tissue anisotropy, (2) incorporating collagen fiber directionality from previous OCT experiments to inform material modeling efforts, and (3) increasing the level of indentation displacement to induce large deformation. Using the recorded experimental force–strain data and computational IFEA methods, regional material parameters for human cervical tissue from a cohort of seven patients with diverse obstetric histories are reported (Table 3).

Material Behavior of Nonpregnant Human Tissue.

Consistent with previous mechanical testing results of human cervical tissue, the indentation force response is time-dependent and nonlinear [15]. The 2D strain field response to indentation (Fig. 6) reveals the mechanical consequence of the collagen fiber architecture. Focusing away from the boundaries, in the midstromal region, directly underneath the indentation location, the strain perpendicular to the fiber direction eyyEXP is larger than the strain parallel to the fibers exxEXP (Fig. 5). A preliminary FEA sensitivity analysis (not shown here) was conducted to understand how the indentation force and strain responses are affected by material parameters, specimen preparation, and boundary conditions. The sensitivity study reveals the magnitude of the separation between exxEXP and eyyEXP, as seen in Fig. 5, is controlled by the dispersion of the collagen fibers. Dispersion of the collagen fibers is characterized by the von-Mises concentration parameter b, where larger values of b indicate more tightly aligned fibers. In computational models of human pregnancy [12], we have found fiber dispersion serves a mechanical function where off-axis fibers resist loading at the internal os.

Overall, the stiffness parameters E and ξ for the equilibrium response are on the order of 1–10 kPa, and ν is 0.22 ± 0.13. The ν indicates a level of compressibility of the material, which is probably due to some level of the interstitial fluid volume decreasing with deformation. Comparing results with previously published uni-axial tension and compression tests, the Young's modulus of the ground substance measured here (E =3.28 kPa) match well with previous work [13,14,16] (E =2 kPa). Hence, spherical indentation is a viable and advantageous method to measure compressive material behavior, where the response is most sensitive to E and ν, without cutting individual uni-axial compression samples and destroying precious samples.

On the other hand, there exists a difference in parameter ξ when comparing material parameters generated from tensile mechanical tests and as reported in Ref. [16]. Fiber stiffness ξ is on the order of 1 kPa measured from indentation and on the order of 10–100 kPa measured from uni-axial tension [13,16]). A slight difference may exist because of small changes to the definition of the fiber strain energy density (Eq. (4)) and the shape of the fiber distribution. The fiber strain energy density is changed to an exponential form from a power-law and the shape of the fiber distribution is changed to a von-Mises from an ellipse. The probable reason for the large difference between measured ξ values is indentation mechanical tests are not sensitive to the material parameter ξ. The fiber stiffness ξ is seen to control the tensile response more than compression response (Fig. 9). Using the ξ value found by indentation, the calculated stress–strain uni-axial material behavior captures the softest specimens tested from previous work [13,16]. Future mechanical investigations will need to couple tensile tests with indentation to fully capture the material behavior for a large tension/compression deformation range.

Heterogeneity of Nonpregnant Human Cervix.

The human cervix has heterogeneous material properties, where the anatomic anterior of the internal os is stiffer than the anatomic anterior of the external os (Figs. 8(a) and 8(b)) and collagen fibers are more tightly aligned in the anterior and posterior quadrants compared to the left and right quadrants (Fig. 8(c)). It is not clear if this pattern of mechanical properties for the cervix is preserved during pregnancy and if a stiffer internal os facilitates the structural loading of the growing fetus. In pregnancy, the internal os is loaded more than the external os [12,17], especially the anterior quadrant. Therefore, it makes sense to find a stiffer internal os compared with the external os in the anterior quadrant. In the continuum material model presented here, the fiber contribution lumps together aligned components. Future studies will tease out the contribution of the individual cellular and ECM components of the cervix.

Each sample in this study was imaged using OCT before mechanical tests to determine the collagen fiber orientation in the plane perpendicular to the inner canal [6]. The material model reported here is informed by the preferential fiber direction of the stroma region measured by OCT, but the dispersion parameter b was left as a fitting parameter determined by the 2D strain field. The magnitude of b measured via indentation mechanical test differed from the value measured via OCT, where the values measured via mechanical test are larger. Further studies are currently being done to understand the correlation between dispersion parameters measured via OCT imaging and mechanical indentation. Regardless of this difference, the trend for fiber dispersion found via indentation is the same as the OCT analysis. The anterior and posterior quadrants of the cervical sample have larger b values, compared with the left and right quadrants, i.e., the anterior and posterior quadrant of the cervix have more tightly aligned fibers.

Experimental and Modeling Error Analysis.

The error in the force–indentation response is a summation of load cell tolerance, the displacement resolution of the material tester, and the buoyancy of the indenter. The error due to buoyancy (Fb<ρgV=2700×π(0.0032×0.003)<0.0001N) is negligible, therefore we left it out of the error calculation.

To analyze the error from the DIC image correlation process, we conducted a rigid body movement study and an enlargement study. In the rigid body study, we used Photoshop (cc 2017, Adobe, Inc., San Jose, CA) to move the selected area randomly to different positions. Using the generated images, we calculated the Lagrange strains using the same setup as we used in our real DIC analyses. The calculated Lagrange strains were less than 10−5, which are several orders less than the experimental indentation strains. For the enlargement study, similarly, we used the Photoshop to enlarge the selected area to specific a ratio (110%, 120%, 130%, and 140%, respectively), and saved them as the new .tif files. We analyzed the generated images using the same DIC setup as we used in our experiments. The errors between the Lagrange strains calculated and the theoretical values are less than 0.05%. Based on this analysis, the errors from the image correlation processes are considered negligible.

For the effects of the glued boundaries, we can clearly see zero displacement on the DIC strain maps for the sample edges (Fig. 4(a)). We accounted for this zero-displacement boundary in the IFEA analysis. Additionally, we have restrained our material model fits to the region directly underneath the indenter tip. Hence, we are confident the glued boundaries did not affect our results.

Optimization Considerations.

Admittedly, although GA has shown advantages in finding the global optima for complex problems [2527], it is still difficult to verify if the search has been trapped into local minima. GA attempts to mitigate this issue by seeding the entire solution space randomly throughout and looking for the global minima. Additionally, by iterating GA for as many as 40 generations, the chance of being trapped into the local optima can be efficiently reduced. In addition, by comparing the exhaustive method for one randomly chosen case, the GA has proven to successfully find the global optima in that case.

Limitations.

Overall, the simplified FEA model of the cervical slices managed to represent the indentation experiments well. However, from the 52 indentation locations on the 13 cervical slices we tested from the 7 patients, the material model fits for 5 indentation locations were not ideal. A representative bad fit (the worst one) is in Fig. 10. The main reasons for bad fitting results are the FEA of the specimens did not capture its complex shape and did not account for the boundary conditions. For these cases, the indentation force response is accurately calculated. However, the strains are not well captured. Although indentation mechanical tests have minimal sample preparation, ill-shaped samples are difficult to model and indentation imposes a complex deformation pattern.

Based on the error field of our validation study (Figs. 7(c,1) and 7(c,2)), normal strains around the indenter are captured well. The material model does a good job of describing the normal strain fields up to the position of 1 mm from the indenter. Although the pattern of shear is captured in the model, the magnitude of shear strain is not well predicted (Fig. 7(c,3)). The shear strain carries a high error away from the indenter. For the area beyond 1 mm from the indenter, the strains are influenced by the irregular boundary conditions, including the shape of the canal and the edge of the glued boundary. Additionally, the chosen fiber model may not accurately reflect the entangled nature of the fiber network. Future work will explore the influence of the boundary and the addition of material parameters that can account for fiber interactions.

Due to the limitation of the current version of the DIC software, we were unable to convert the DIC coordinate frame to a circular coordinate frame. In future studies, we will code our own method to enable visualization of strain fields in an appropriate coordinate frame for our biologic samples.

Finally, due to the scarcity of samples, this study represents a small cohort of seven women. Hence, this study is not statistically powered to make conclusions on the influence of patient demographics on the mechanical behavior of the cervix. Further mechanical tests using this methodology are being conducted on a larger cohort of specimens to understand the effects of parity, age, race, and obstetric history.

Conclusions

The mechanical properties of the nonpregnant human cervix are anisotropic and nonlinear, as characterized by spherical indentation mechanical tests coupled with strain field measurements. A fiber network composite material model captures salient features of the equilibrium material response of the tissue, where experimental conditions are most sensitive to the material parameters associated with the ground substance [E, ν] and the measure of fiber dispersion [b]. Based on comparing material properties between indentation locations within a cervical sample, the cervix is a heterogeneous organ, where the anterior quadrant of its internal os is stiffer than the anterior quadrant of its external os. When comparing anatomical quadrants, the tissue in the anterior and posterior quadrants tends to be more anisotropic than tissue from the left and right quadrants. The reason for this pattern of mechanical properties remains to be determined. Finally, there is no significant difference found between the tissue samples from women with lower parity and the tissue from women with higher parity. In conclusion, indentation modality is a viable and nondestructive method to mechanically interrogate the cervix, and it offers a potential way to measure the in vivo mechanical properties of the cervix during pregnancy in hopes of assessing the risk of preterm birth.

Acknowledgment

We acknowledge computing resources from Columbia University's Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR) Contract C090171, both awarded April 15, 2010.

We want to pay special acknowledgment to the pioneering work of Dr. Y. C. Fung, for whom this work is dedicated to. The reproductive biomechanics work featured in this study is directly influenced by his work in nonlinear elasticity, viscoelasticity, and other continuum mechanical models applied to living tissues. We aspire to solve the problem of preterm birth one day by the careful and meaningful application of biomechanics to the study of human pregnancy. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Funding Data

  • National Institutes of Health (Grant No. R01HD091153 to KM; Funder ID: 10.13039/100000002).

  • NIH Research Facility Improvement (Grant No. 1G20RR030893-01; Funder ID: 10.13039/501100000272).

  • New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR) (Funder ID: 10.13039/100011036).

Appendix: Correlation of Finite Element Analysis and DIC Coordinates

The 2D strain fields in this work were measured in the coordinate frame of the DIC cameras (Fig. 4(a)), and the strains were denoted as exxDIC,eyyDIC, and exyDIC. Whenever the experimental strains were compared against strains calculated with the FE model, either for fitting or validation purposes, the experimental strain data were rotated using the following equations: 
{exxDIC,cor=exxDICcos2θ+eyyDICsin2θ+exyDICsinθcosθeyyDIC,cor=exxDICsin2θ+eyyDICcos2θ+exyDICsinθcosθexyDIC,cor=2(eyyDICexxDIC)sinθcosθ+exyDIC(cos2θsin2θ)
(A1)
where θ was the angle difference between FEA and DIC coordinates (Fig. 7(a)), exxDIC,cor and eyyDIC,cor were the Lagrange strains in the transformed coordinate frame (FE). For convenience, we omitted the superscript, “cor” and express “DIC” as “EXP” for the experimental data. In other words, we used exxEXP, eyyEXP, and exyEXP to represent the experimental strains in the FE coordinate frame, as done in the optimization routine in Sec. 2.6 and validation in Sec. 2.7.

In this way, the strains obtained from the DIC results and the FEA simulations could be compared. exx and eyy are directly the Lagrangian strains in the circumferential and radial direction, respectively, under the indenter. Then, only strains at the four equilibrium levels (exxEXPi, eyyEXPi, i=1,2,3,4) and the corresponding force response data (FEXPi, i=1,2,3,4) were selected for the data fitting process.

References

References
1.
Myers
,
K. M.
,
Feltovich
,
H.
,
Mazza
,
E.
,
Vink
,
J.
,
Bajka
,
M.
,
Wapner
,
R. J.
,
Hall
,
T. J.
, and
House
,
M.
,
2015
, “
The Mechanical Role of the Cervix in Pregnancy
,”
J. Biomech.
,
48
(
9
), pp.
1511
1523
.
2.
Feltovich
,
H.
,
Hall
,
T.
, and
Berghella
,
V.
,
2012
, “
Beyond Cervical Length: Emerging Technologies for Assessing the Pregnant Cervix
,”
Am. J. Obstet. Gynecol.
,
207
(
5
), pp.
345
354
.
3.
Vink
,
J.
, and
Feltovich
,
H.
,
2016
, “
Cervical Etiology of Spontaneous Preterm Birth
,”
Semin. Fetal Neonat. Med.
,
21
(
2
), pp.
106
112
.
4.
WHO
,
2018
, “
World Health Organization Fact Sheet No. 363
,” World Health Organization, accessed Dec. 19, https://www.who.int/news-room/fact-sheets/detail/preterm-birth
5.
Bauer
,
M.
,
Mazza
,
E.
,
Nava
,
A.
,
Zeck
,
W.
,
Eder
,
M.
,
Bajka
,
M.
,
Cacho
,
F.
,
Lang
,
U.
, and
Holzapfel
,
G. A.
,
2007
, “
In Vivo Characterization of the Mechanics of Human Uterine Cervices
,”
Ann. N. Y. Acad. Sci.
,
1101
(
1
), pp.
186
202
.
6.
Yao
,
W.
,
Gan
,
Y.
,
Myers
,
K. M.
,
Vink
,
J. Y.
,
Wapner
,
R. J.
, and
Hendon
,
C. P.
,
2016
, “
Collagen Fiber Orientation and Dispersion in the Upper Cervix of Non-Pregnant and Pregnant Women
,”
PLoS One
,
11
(
11
), p.
e0166709
.
7.
Zork
,
N. M.
,
Myers
,
K. M.
,
Yoshida
,
K.
,
Cremers
,
S.
,
Jiang
,
H.
,
Ananth
,
C. V.
,
Wapner
,
R. J.
,
Kitajewski
,
J.
, and
Vink
,
J.
,
2014
, “
A Systematic Evaluation of Collagen Cross-Links in the Human Cervix
,”
Am. J. Obstet. Gynecol.
,
212
(
3
), pp.
1
8
.
8.
House
,
M.
,
Kaplan
,
D.
, and
Socrate
,
S.
,
2009
, “
Relationships Between Mechanical Properties and Extracellular Matrix Constituents of the Cervical Stroma During Pregnancy
,”
Semin. Perinatol.
,
33
(
5
), pp.
300
307
.
9.
Vink
,
J. Y.
,
Qin
,
S.
,
Brock
,
C. O.
,
Zork
,
N. M.
,
Feltovich
,
H. M.
,
Chen
,
X.
,
Urie
,
P.
,
Myers
,
K. M.
,
Hall
,
T. J.
,
Wapner
,
R.
,
Kitajewski
,
J. K.
,
Shawber
,
C. J.
, and
Gallos
,
G.
,
2016
, “
A New Paradigm for the Role of Smooth Muscle Cells in the Human Cervix
,”
Am. J. Obstet. Gynecol.
,
215
(
4
), pp.
478.e1
478.e11
.
10.
Aspden
,
R.
,
1988
, “
Collagen Organization in the Cervix and Its Relation to Mechanical Function
,”
Collagen Relat. Res.
,
8
(
2
), pp.
103
112
.
11.
Weiss
,
S.
,
Jaermann
,
T.
,
Schmid
,
P.
,
Staempfli
,
P.
,
Boesiger
,
P.
,
Niederer
,
P.
,
Caduff
,
R.
, and
Bajka
,
M.
,
2006
, “
Three-Dimensional Fiber Architecture of the Nonpregnant Human Uterus Determined Ex Vivo Using Magnetic Resonance Diffusion Tensor Imaging
,”
Anat. Rec., Part A
,
288
(
1
), pp.
84
90
.
12.
Fernandez
,
M.
,
House
,
M.
,
Jambawalikar
,
S.
,
Zork
,
N.
,
Vink
,
J.
,
Wapner
,
R.
, and
Myers
,
K.
,
2015
, “
Investigating the Mechanical Function of the Cervix During Pregnancy Using Finite Element Models Derived From High-Resolution 3D MRI
,”
Comput. Methods Biomech. Biomed. Eng.
,
19
(
4
), pp.
404
417
.
13.
Myers
,
K.
,
Socrate
,
S.
,
Paskaleva
,
A.
, and
House
,
M.
,
2010
, “
A Study of the Anisotropy and Tension/Compression Behavior of Human Cervical Tissue
,”
ASME J. Biomech. Eng.
,
132
(
2
), p.
021003
.
14.
Myers
,
K.
,
Paskaleva
,
A.
,
House
,
M.
, and
Socrate
,
S.
,
2008
, “
Mechanical and Biochemical Properties of Human Cervical Tissue
,”
Acta Biomater.
,
4
(
1
), pp.
104
116
.
15.
Yao
,
W.
,
Yoshida
,
K.
,
Fernandez
,
M.
,
Vink
,
J.
,
Wapner
,
R.
,
Ananth
,
C.
,
Oyen
,
M.
, and
Myers
,
K.
,
2014
, “
Measuring the Compressive Viscoelastic Mechanical Properties of Human Cervical Tissue Using Indentation
,”
J. Mech. Behav. Biomed. Mater.
,
34
, pp.
18
26
.
16.
Myers
,
K. M.
,
Hendon
,
C. P.
,
Gan
,
Y.
,
Yao
,
W.
,
Yoshida
,
K.
,
Fernandez
,
M.
,
Vink
,
J.
, and
Wapner
,
R. J.
,
2015
, “
A Continuous Fiber Distribution Material Model for Human Cervical Tissue
,”
J Biomech
,
48
(
9
), pp.
1533
1540
.
17.
Westervelt
,
A. R.
,
Fernandez
,
M.
,
House
,
M.
,
Vink
,
J.
,
Nhan-Chang
,
C.-L.
,
Wapner
,
R.
, and
Myers
,
K. M.
,
2017
, “
A Parameterized Ultrasound-Based Finite Element Analysis of the Mechanical Environment of Pregnancy
,”
ASME J. Biomech. Eng.
,
139
(
5
), p.
051004
.
18.
Feng
,
Y.
,
Okamoto
,
R. J.
,
Namani
,
R.
,
Genin
,
G. M.
, and
Bayly
,
P. V.
,
2013
, “
Measurements of Mechanical Anisotropy in Brain Tissue and Implications for Transversely Isotropic Material Models of White Matter
,”
J. Mech. Behav. Biomed. Mater.
,
23
, pp.
117
132
.
19.
Zhang
,
M.
,
Zheng
,
Y.
, and
Mak
,
A. F.
,
1997
, “
Estimating the Effective Young's Modulus of Soft Tissues From Indentation Testsnonlinear Finite Element Analysis of Effects of Friction and Large Deformation
,”
Med. Eng. Phys.
,
19
(
6
), pp.
512
517
.
20.
Namani
,
R.
,
2012
, “
Elastic Characterization of Transversely Isotropic Soft Materials by Dynamic Shear and Asymmetric Indentation
,”
ASME J. Biomech. Eng.
,
134
(
6
), p.
061004
.
21.
Elkin
,
B. S.
, and
Morrison
,
B.
,
2013
, “
Viscoelastic Properties of the P17 and Adult Rat Brain From Indentation in the Coronal Plane
,”
ASME J. Biomech. Eng.
,
135
(
11
), p.
114507
.
22.
Lake
,
S. P.
, and
Barocas
,
V. H.
,
2012
, “
Mechanics and Kinematics of Soft Tissue Under Indentation Are Determined by the Degree of Initial Collagen Fiber Alignment
,”
J. Mecha. Behav. Biomed. Mater.
,
13
, pp.
25
35
.
23.
Han
,
L.
,
Noble
,
J. A.
, and
Burcher
,
M.
,
2003
, “
A Novel Ultrasound Indentation System for Measuring Biomechanical Properties of In Vivo Soft Tissue
,”
Ultrasound Med. Biol.
,
29
(
6
), pp.
813
823
.
24.
Holland
,
J. H.
,
1992
,
Adaptation in Natural and Artificial Systems: An Introductory Analysis With Applications to Biology, Control, and Artificial Intelligence
,
MIT Press
,
Cambridge, MA
.
25.
Saga
,
M.
, and
Vasko
,
M.
,
2007
, “
Solution of Mechanical Systems With Uncertainty Parameters Using IFEA
,”
Communications
,
11
(
2
), pp.
19
27
.http://www.utc.sk/komunikacie/archiv/2009/2/2_2009en.pdf
26.
Chawla
,
A.
,
Mukherjee
,
S.
, and
Karthikeyan
,
B.
,
2009
, “
Characterization of Human Passive Muscles for Impact Loads Using Genetic Algorithm and Inverse Finite Element Methods
,”
Biomech. Model. Mechanobiol.
,
8
(
1
), pp.
67
76
.
27.
He
,
Y.
,
Guo
,
D.
, and
Chu
,
F.
,
2001
, “
Using Genetic Algorithms and Finite Element Methods to Detect Shaft Crack for Rotor-Bearing System
,”
Math. Comput. Simul.
,
57
(
1–2
), pp.
95
108
.
28.
Tadepalli
,
S. C.
,
Erdemir
,
A.
, and
Cavanagh
,
P. R.
,
2011
, “
Comparison of Hexahedral and Tetrahedral Elements in Finite Element Analysis of the Foot and Footwear
,”
J. Biomech.
,
44
(
12
), pp.
2337
2343
.
29.
Myers
,
K.
, and
Ateshian
,
G.
,
2014
, “
Interstitial Growth and Remodeling of Biological Tissues: Tissue Composition as State Variables
,”
J. Mech. Behav. Biomed. Mater.
,
29
, pp.
544
556
.
30.
Gasser
,
T.
,
Ogden
,
R.
, and
Holzapfel
,
G.
,
2006
, “
Hyperelastic Modelling of Arterial Layers With Distributed Collagen Fibre Orientations
,”
J. R. Soc. Interface
,
3
(
6
), pp.
15
35
.