Abstract

Vaginal childbirth is the final phase of pregnancy when one or more fetuses pass through the birth canal from the uterus, and it is a biomechanical process. The uterine active contraction, causing the pushing force on the fetus, plays a vital role in regulating the fetus delivery process. In this project, the active contraction behaviors of muscle tissue were first modeled and investigated. After that, a finite element method (FEM) model to simulate the uterine cyclic active contraction and delivery of a fetus was developed in ls-dyna. The active contraction was driven through contractile fibers modeled as one-dimensional truss elements, with the Hill material model governing their response. Fibers were assembled in the longitudinal, circumferential, and normal (transverse) directions to correspond to tissue microstructure, and they were divided into seven regions to represent the strong anisotropy of the fiber distribution and activity within the uterus. The passive portion of the uterine tissue was modeled with a Neo Hookean hyperelastic material model. Three active contraction cycles were modeled. The cyclic uterine active contraction behaviors were analyzed. Finally, the fetus delivery through the uterus was simulated. The model of the uterine active contraction presented in this paper modeled the contractile fibers in three-dimensions, considered the anisotropy of the fiber distribution, provided the uterine cyclic active contraction and propagation of the contraction waves, performed a large deformation, and caused the pushing effect on the fetus. This model will be combined with a model of pelvic structures so that a complete system simulating the second stage of the delivery process of a fetus can be established.

1 Introduction

Normal biomechanical function and regulation of the uterus are vital for the childbirth process. Poor quality uterine contractions contribute to many labor problems, including labor dystocia (a delay in labor progression), which has proven to be a significant reason for the rising Cesarean section rate [1]. In addition, the interaction of the uterus, the fetus, and the maternal pelvis results in both the normal cardinal movements of labor and abnormalities of descent that can result in a shoulder dystocia (delay in delivering the shoulders after the head delivers). Therefore, it is important to be able to model uterine contraction behavior during the delivery process in order to better understand the pathomechanics of the process itself and its relationship to injuries that can occur in both the mother and the infant.

The uterine wall consists of three layers: the endometrium, the myometrium, and the perimetrium (from inside to outside). The myometrium, the thickest layer, is composed of millions of smooth muscle cells (SMC). Such SMCs are organized into fiber-like constructs with different orientations, and this provides the complex contraction behavior of the uterus [2]. The longitudinal and circumferential direction are the two main fiber directions inside of the uterine wall [35]. Recent studies have shown that the fibers inside of the uterine wall for both nonpregnant and pregnant women are not perfectly parallel to the uterine surface but are inclined within the uterine wall, which results in a three-dimensional organization [6,7]. The distribution of the fibers shows strong anisotropy and heterogeneity in different regions of the uterus, and such variations have a significant effect on uterine mechanical behaviors [8,9]. The contraction behaviors of the uterus are summarized in Grimm's review [10]. The typical contraction pattern of the uterus resembles a bell-like curve [11], which takes on average 60 s for one cycle [12]. The active contraction is a force wave that is initiated from the fundus and then propagates to the lower part of the uterus, which takes 10–20 s [13]. While the contraction in different regions starts at various times, the peak force occurs almost at the same time throughout the uterus, and then relaxation occurs synchronously [13].

Due to ethical issues, human experiments during labor and delivery are difficult to conduct. Computational modeling has become a promising way to investigate the childbirth process. For example, Myers have done several works about characterization of mechanical properties of maternal tissues [1416], constitutive modeling [17,18], and finite element analysis of prebirth process [19,20]. Oyen developed FEM models to investigate the rupture of the uterus [21] and fibrous tissue during labor [22,23]. There are also some other works about simulating the movement of the fetus through the pelvis structures without uterus [24,25]. However, simulation of uterine active contraction is still a relatively new research area. In 2009 and 2017, Buttin et al. [26,27] established an FEM model of the female reproductive system containing the pelvis, uterus, and fetus to investigate how the different forces influence the fetus delivery process. However, in these simulations, the uterus's active contraction was not simulated directly, but rather an external force was applied through the uterus to the rigid fetal model to simulate what was assumed to be the resulting force. In 2015, Cochran et al. [28] developed a three-dimensional electromechanical coupling model of the uterus to calculate the mechanical and electrical properties and the peak uterine pressure. A similar model was developed by Yochum et al. [29] in 2016, where the mechanical and electrical behaviors were analyzed on multiple scales. However, the uterine geometries in these simulation models were highly simplified due to their complexities, and the calculated strain was too small to represent the large deformation of the uterus during delivery. In 2019, Pouca, et al. [30] developed a chemo-mechanical coupling model to simulate uterine contraction and fetal delivery. However, the fibers only had two directions, the longitudinal and circumferential direction, and the effects of anisotropy of fibers distribution were not considered. Overall, the existing models in this area still have many limitations.

In this project, the active contraction behaviors of a muscle tissue were first modeled and investigated. Based on this, a FEM uterus model was developed in ls-dyna to simulate the uterine active contraction, where the contraction forces were driven by contractile fibers. Multiple mechanical contraction properties of the uterus were included in this model. Following the development of the model, it was used to demonstrate its utility in simulating pushing to deliver a fetus.

2 Materials and Methods

2.1 Muscle Tissue.

A generalized structure of a unit of muscle tissue was created within ls-dyna, and its mechanical behaviors was investigated as the first step in this project (Fig. 1). The passive portion of the muscle—representing the noncontractile components of the cells and extracellular connective structures—was modeled as a 100 × 10 × 10 mm solid (red portion in Fig. 1(a)) for convenience in analysis. The active, contractile fibers of the muscle were modeled using four truss elements, each with a length of 100 mm and a cross-sectional area of 100 mm2 (blue components in Fig. 1(b)). Each of the four contractile fibers was placed along one of the four long edges of the passive solid. The passive solid portion and the four active beams were coupled by sharing the nodes between the truss elements and the solid elements. There were no loading or boundary conditions applied to the muscle tissue model.

Fig. 1
Geometrical model the muscle tissue (a) and activation level for contractile fibers (b)
Fig. 1
Geometrical model the muscle tissue (a) and activation level for contractile fibers (b)
Close modal

The Hill material [31] was used for the beams, while an elastic material model (E = 2 kPa, v = 0.49) was initially used for the passive solid. The activation level, with a range of 0–1, for the Hill model was assumed as shown in Fig. 1(b). The normalized force–length curve (SVS) and force–velocity curve (SVR), which describe fundamental responses of muscle and have previously been shown to be appropriate to describe smooth muscle's behavior [32,33], were obtained from the literature [34]. The maximum isometric stress (PIS) was initially set at 50 Pa, which was a trial value to generate a relatively small deformation based on the Young's modulus of the muscle tissue, and the damping coefficient (DMP) was 0.001, which was close to 0.004 used previously in the literature for muscle modeling [34]. Such a set of parameters was not selected to represent a specific muscle tissue. There were two purposes to modeling such a generalized muscle tissue. The first was to verify whether the bonding of the active contractile fibers—using the Hill model—and a passive portion of the solid elements would work to simulate a muscle's active contraction and relaxation behaviors. The second was to investigate and understand how each of the parameters affects the muscle tissue's mechanical behaviors (such as the contraction velocity, and contraction capacity) quantitively and qualitatively. Therefore, a parametric analysis was also conducted. The range for the parameters is provided in Table 1. The Young's modulus values selected of 2, 20, and 40 kPa were within the normal range of values for muscle tissue. However, the range of values for DMP and PIS for muscle tissues, especially smooth muscle tissues, were unavailable from the literature. As such, a range was assumed as shown in Table 1.

Table 1

Range of parameters investigated for generalized muscle tissue

ParameterRange
Young's modulus (E)2, 20, 40 kPa
Isometric maximum stress (PIS)5 × 10−5, 5 × 10−4, 5 × 10−3, 5 × 10−2 MPa
Damping coefficient (DMP)1 × 10−4, 1 × 10−3, 1 × 10−2
ParameterRange
Young's modulus (E)2, 20, 40 kPa
Isometric maximum stress (PIS)5 × 10−5, 5 × 10−4, 5 × 10−3, 5 × 10−2 MPa
Damping coefficient (DMP)1 × 10−4, 1 × 10−3, 1 × 10−2

2.2 Uterine Model.

The geometric models of the uterus and the fetus (Fig. 2) were designed in NX software based on real-time magnetic resonance imaging of the uterus and fetus during the second stage of labor [35]. The inner diameter of the cervix was 100 mm, and the thickness of the uterine wall was set to be 10 mm (Fig. 2(a)). The length and maximum width of the uterus were 349 mm and 156 mm, respectively (Fig. 2(b)). The fetus model included three components: fetal head, neck, and body. The fetal head was a sphere with a diameter of 80 mm, and the length of the whole fetus (with arms and legs tucked against the body) was about 330 mm. All dimensions were representative of values that might be observed at 40 weeks of gestation.

Fig. 2
Geometrical models of the fetus and the uterus
Fig. 2
Geometrical models of the fetus and the uterus
Close modal

Meshwork of the uterus (Fig. 3(a)) and fetus (Fig. 3(b)) was established in Hypermesh software. Hexahedral elements and tetrahedral elements were used for the uterus and the fetus, respectively. The hexahedral elements were chosen and one layer of such elements was defined in the normal direction for the uterine wall because it is more convenient to bond with the contractile fibers. There were 3636 nodes and 1795 elements for the uterus structure, and 7367 nodes and 33,034 elements for the fetus in total. The mesh for these two structures was then exported to ls-dyna to assign material properties and boundary conditions and to perform the analysis. The contact between the uterus and the fetus was defined as a surface-to-surface contact.

Fig. 3
Mesh of the uterus (a) and the fetus (b)
Fig. 3
Mesh of the uterus (a) and the fetus (b)
Close modal

The fibers inside of the uterine wall have distinct directional organization based on experimental findings [7]; therefore, it was important to model the contractile fibers in three-dimensions in a FEM simulation model, not only in two-dimensions—with fibers in the longitudinal and circumferential directions parallel to the surface of the uterine muscle tissue—as has been reported in some previous studies [30]. The three directions of fibers can be used to investigate the effects of fiber orientation on uterine mechanical behaviors more accurately. From the mechanical perspective, any direction of force can be divided into three basic forces along the x, y, and z-coordinate axes. Therefore, the simulation model in this study created the contractile fibers inside of the uterine muscle wall in three directions–the longitudinal, circumferential, and normal directions as shown in Figs. 4(a)4(c).

Fig. 4
Contractile fibers inside of the uterine wall: longitudinal (a), circumferential (b), normal (c), and the combination of these three directions of fibers (d)
Fig. 4
Contractile fibers inside of the uterine wall: longitudinal (a), circumferential (b), normal (c), and the combination of these three directions of fibers (d)
Close modal

In addition, the distribution for both the orientation and content of the fibers changes dramatically in different regions of the uterus [5,8,36]. Also, the contraction of the uterus has been shown to be a force wave starting from the fundus region and propagating to the lower segments of the uterus, with a higher contraction intensity in the fundal region [13]. To consider such anisotropies and heterogeneities of the fibers, the three types of contractile fibers in the longitudinal, circumferential, and normal directions modeled in the study were all divided into seven regions (Fig. 4) so that the uterine contraction behaviors can be better modeled. These seven regions were independent of one another. The seven independent regions result in a model in which it is possible to adjust the fiber distribution by changing the cross-sectional area in different regions individually while not influencing the fibers in other regions. In addition, in each region, the three directions of fibers were also independent of each other, which means that adjusting the cross-sectional area for one direction of fiber will not affect the fibers in other directions in the same region. This will impact the resultant direction of force generated by the contraction of each muscle element. The model thus has great flexibility to adjust the distribution of the three orientations of fibers within the different regions. The passive portion of the uterine tissue was constant throughout the model. Figure 4 shows the seven regions of the longitudinal (Fig. 4(a)), circumferential (Fig. 4(b)), and normal direction fibers (Fig. 4(c)) along with the combination of these three types of fibers (Fig. 4(d)). The seven different colors for each direction of fiber (Figs. 4(a)4(c)), were used to differentiate the fibers of the seven different regions. The seven regions were arranged from the top to the bottom of the uterus. At the interfaces between every two regions, they share the same element nodes.

In this uterine simulation model, the average size of a solid element was 10 mm × 10 mm × 10 mm (Fig. 3(a))—with a cross-sectional area of 100 mm2. For each solid element, there are twelve edges that are all bonded with fibers (Fig. 4(d)): four in the longitudinal direction (Fig. 4(a)); four in the circumferential direction (Fig. 4(b)); and four in the normal direction (Fig. 4(c)). The cross-sectional area for each longitudinal fiber was 20 mm2 (Fig. 5), while each longitudinally oriented edge of a solid element was shared by two elements (Fig. 3(a)). As a result, the contribution of each longitudinal fiber for a single solid element would be 10 mm2. The total cross section area from these four associated longitudinal fibers is thus 40 mm2; therefore, the content of longitudinal fibers would be approximately 40% for a random solid element (Fig. 5). The reason for the approximation is that the complex shape of the uterus results in some variation in the element shape from a perfect cube. Similarly, the contribution of the four 20 mm2 circumferential fibers to each element was around 40% of the cross-sectional area in the fundus region. The cross-sectional area of each normal fiber was set as 6 mm2 in the model, in order to represent the fact that the fibers within the uterine wall are not perfectly parallel to the uterine wall surface but are inclined with small angles [7] Each edge on which a normal fiber was located is shared with four surrounding solid elements (Fig. 3(a)), so the contribution of the four normal fibers was 6% of the cross-sectional area of the element. As a result, the volume content of all of the fibers in total was about 86%, which is close to the experimental finding that—in a fundal region of a pregnant uterus' myometrium—smooth muscle cells that equate to the fibers account for approximately 80% of the total tissue volume [7]. Also, in the midcorpus of the uterus, the longitudinal direction fibers were the dominant ones [8]. In this model, the cross-sectional area of the circumferential fibers was set as 6 mm2 in the midcorpus, with the other two direction fibers unchanged. As a result, the total fiber content in the midcorpus was 58%, which was also comparable to the range of fiber content measured in the midcorpus for the pregnant and unpregnant uterus of 40% and 63%, respectively [37].

Fig. 5
A random element showing the muscle tissue volume and four fibers in the longitudinal direction. L, C, and N represent the longitudinal, circumferential, and normal directions, respectively.
Fig. 5
A random element showing the muscle tissue volume and four fibers in the longitudinal direction. L, C, and N represent the longitudinal, circumferential, and normal directions, respectively.
Close modal

The contractile fibers were modeled as truss elements with the Hill material model [31]. The SVS and SVR were kept the same as the functions used in modeling the general smooth muscle tissue [34]. A PIS of 1.0 MPa and a DMP of 100 were used for the model of uterine muscle. The parameters of PIS and DMP for the Hill model were assumed based on the degree of tissue contraction that is needed to push the fetus from the uterus with an overall delivery time that aligns with clinical observation. For PIS, which again is the maximum isometric stress that can be developed by the tissue, the selected value of 1 MPa is comparable to values of 0.5–1 MPa used in prior muscle modeling described in the literature [38]. Clinically, delivery time can vary significantly, from 2 min to 200 min [39]. In this simulation, it was assumed that the fetus would be delivered within 3 min. The activation-level curves for the seven regions of contractile fibers were assumed as shown in Fig. 6. The starting time for the activation of contraction for region one (fundal region) of the uterus was at time 0 s. Each of the subsequent regions began contracting 2 s later than the previous, higher region to simulate the propagation of the contraction wave. The propagation of a contraction wave initiated from the fundal region to the middle and then lowest region of the uterus takes about 10–20 s to finish, according to experimental findings [13] and has been modeled as 12 s in previous simulations [28,29]. In this model, the uterus was divided into seven regions, and the delay for each region was 2 s so that the time period for the propagation was 12 s in total. All seven regions within the uterus reached a fully activated status and then began to relax synchronously, which was based on the physiologic behavior of the uterus [13]. Three uterine active contraction cycles were modeled in this study (Fig. 6).

Fig. 6
Activation level curves for the seven regions of contractile fibers
Fig. 6
Activation level curves for the seven regions of contractile fibers
Close modal

For the material properties of the uterus and the fetus, the passive portion of the uterus muscle, representing the noncontractile portions of the tissue, was modeled as a hyperelastic material with a Neo Hookean model (C1 = 0.03 MPa) [27]. The fetus was modeled as a rigid body. The solid hexahedral elements for the uterus were coupled with the truss elements for the contractile fibers by sharing the nodes (Fig. S1 available in the Supplemental Materials on the ASME Digital Collection). The cervix was fixed in the space as the boundary condition.

3 Results

3.1 Active Contraction Behaviors of the Generalized Muscle Tissue.

The contraction and relaxation of the muscle tissue along the x-axis are shown in Fig. 7. The muscle started to contract from the two ends toward the middle once the model was activated. At t = 2 s, the muscle model had shortened by 5.74 mm (strain = 5.74%) (Fig. 7(b)). At t = 40 s, the solid reached its contraction capacity with a strain of 8.09% (Fig. 7(c)). When the activation level began to decrease (t = 40 s), the muscle started to relax and return to its original state. At t = 100 s, the muscle tissue had completely returned to the initial shape (Fig. 7(d)). The mechanical behaviors of the muscle tissue were regulated by both the contractile fibers and the mechanical properties of the passive portion of the muscle tissue. The active contraction force generated by the contractile fibers made the muscle tissue contract, while the passive portion of the muscle tissue represented by the solid elements resisted such contraction and supported the complete relaxation of the muscle tissue unit to its original length.

Fig. 7
Deformation for the muscle tissue at times of 0 s (a), 2 s (b), 40 s (c), and 100 s (d)
Fig. 7
Deformation for the muscle tissue at times of 0 s (a), 2 s (b), 40 s (c), and 100 s (d)
Close modal

The contraction capacity (the maximum contraction deformation that a muscle tissue unit can develop) and contraction speed are two of the most important variables to describe a muscle tissue's mechanical behaviors. The parameter PIS (maximum isometric stress) within the Hill model was found to significantly influence the contraction capacity. Specifically, for all three Young's modulus values used to characterize the passive portion of the muscle (2, 20, and 40 kPa), the maximum strain produced within the solid increased with an increase in PIS (Fig. 8(a)). As would be expected, lower modulus values allowed more deformation during contraction for each level of PIS (Fig. 8(a))—the maximum isometric stress that can be generated will cause greater deformation for tissue with a lower modulus.

Fig. 8
Relationship between the maximum strain of a muscle tissue and PIS (a) and curves describing the deformation due to contraction as a function of time with different values of DMP (b)
Fig. 8
Relationship between the maximum strain of a muscle tissue and PIS (a) and curves describing the deformation due to contraction as a function of time with different values of DMP (b)
Close modal

Likewise, the contraction speed of a muscle can be significantly affected by the DMP coefficient in the Hill material model. The contraction speed, which was the slope of the curves in Fig. 8(b), increased with a decrease in DMP. DMP represents a viscous component of the tissue, and greater viscosity will result in a longer time needed for the muscle unit to contract. The maximum amount of displacement was still governed by the combination of Young's modulus and PIS—which was demonstrated by the fact that the curves in Fig. 8(b) for damping levels of 0.001 and 0.0001 both plateaued at the same level. As the rate of contraction was much lower at the damping level of 0.01, the amount of contraction that was achieved was less than the maximum possible and was limited by the time of activation.

3.2 Active Contraction Behaviors of the Uterus.

3.2.1 Propagation of the Contraction Wave.

The simulation results for the propagation of the contraction wave through the uterus are shown in Fig. 9. Each row of images in Fig. 9 represents a single contraction cycle in the model. For the first cycle (first row of Fig. 9), there was no deformation at t = 0 s. After that, the uterus started to contract from the top of the fundus region and the contraction wave started to propagate to the lower portion of the uterus. By t = 2 s, the whole fundus region of the uterus had contracted. By t = 8 s, the first four regions had started to contract. By t =12 s, all of the seven regions within the uterus had contracted. Once the uterus reached the maximum deformation, the uterus started to relax due to the decrease in the activation level. The initial strain for the second and third contraction cycles was not zero because there was some deformation remaining from the previous cycle—the passive component of the tissue had not completely returned to its original state. Such a phenomenon occurred because the damping effect on the uterus requires a longer time for the uterus to recover than allowed due to the beginning of the next contraction cycle. In each of the following cycles, a new contraction wave propagated to the lower part of the uterus. The overall strain within the uterus increased with each subsequent contraction cycle. The model successfully simulated the propagation phenomenon of the contraction wave within the uterus.

Fig. 9
Propagation of the contraction wave for three contraction cycles showing the strain within the uterine tissue
Fig. 9
Propagation of the contraction wave for three contraction cycles showing the strain within the uterine tissue
Close modal

3.2.2 The Element Stress and Stress Distribution Within the Uterus.

The stress of an element, picked in the fundus region of the uterus, and distribution of the stress within the uterus through the three contraction cycles are shown in Fig. 10. As the precise amount of stress cannot be validated against experimental data, the qualitative pattern of stress was the key result. The element stress increased significantly due to the rise of the activation level and decreased slightly due to the decline of the activation level within each cycle (Fig. 10(a)). The element stress continued to increase with the number of contraction cycles (Fig. 10(a)). For the distribution of the stress within the uterus, the stress was almost evenly distributed for the middle and lower part of the uterus, while the largest stress occurred in the fundus region (Fig. 10(b)).

Fig. 10
Stress for an element in the fundus region (a) and stress nephogram for the uterus with a scale bar of 0–0.025 MPa (b). The three red points in (a) represent the three time points for the stress nephogram for the uterus for the first cycle (first row in (b)). The three blue and three green points were the time points picked for the second and third cycle, respectively.
Fig. 10
Stress for an element in the fundus region (a) and stress nephogram for the uterus with a scale bar of 0–0.025 MPa (b). The three red points in (a) represent the three time points for the stress nephogram for the uterus for the first cycle (first row in (b)). The three blue and three green points were the time points picked for the second and third cycle, respectively.
Close modal

The mechanical response of elements in different positions of the uterus during the simulation were investigated. Three elements, in the fundus, middle, and lower position of the uterus (Fig. S2 available in the Supplemental Materials on the ASME Digital Collection), were picked. The element in the fundus region had the largest displacement, while the elements in the middle and lower positions had the second and smallest displacement during the contraction, respectively (Fig. 11(a)). The element in the fundus region had much larger stress and strain compared to the elements in the middle and lower positions (Figs. 11(b) and 11(c)). The middle element had slightly larger stress and strain compared to the element in the lower position (Figs. 11(b) and 11(c)). In addition, there were slight phase delays for the three curves in each part of Fig. 11 because the contraction happened first in the top region of the uterus and then propagated to the middle and the lower parts.

Fig. 11
Curves of displacement (a), stress (b), and strain (c) for elements in fundus, middle, and lower positions of the uterus
Fig. 11
Curves of displacement (a), stress (b), and strain (c) for elements in fundus, middle, and lower positions of the uterus
Close modal

3.2.3 Parametric Analysis.

As in the previous assessment of the generic muscle contracting unit, the PIS and DMP values in the Hill material model were found to significantly affect the contraction behaviors of the uterus. The results for the element in the fundus region were taken as an example. For the same PIS, the contraction velocity (the slope of the displacement versus time curve) decreased with an increase of the DMP value (Fig. 12(a)). For the same DMP value, the largest displacement of the element in the fundus (representing the contraction capacity or largest deformation of the uterus), the contraction velocity, and the element stress and strain all increased significantly with an increase of PIS (Fig. 12). The larger DMP was found to cause a longer time delay for the elements to reach their largest contraction responses for the same PIS. This can be demonstrated by the fact that the displacement, along with the element stress and strain, reached their maximum values quickly in the first contraction cycle for the smallest DMP, while they continued to increase in the second and third contraction cycles for larger DMP cases (Fig. 12). The elements in the middle and lower positions had similar results (Fig. S3 available in the Supplemental Materials).

Fig. 12
Effects of parameters (PIS & DMP) on element's displacement (a), stress (b), and strain (c) for the element in the fundus location
Fig. 12
Effects of parameters (PIS & DMP) on element's displacement (a), stress (b), and strain (c) for the element in the fundus location
Close modal

In conclusion, the mechanical responses of the elements varied in different positions within the uterus. The elements in the fundus region had larger displacements along with larger stress and strain values than elements in the lower part of the uterus. There was a delay for such responses in the lower elements due to the propagation of the contraction wave that was initiated from the top region (fundus) of the uterus. In addition, the contraction capacity or the largest deformation that a uterus can develop was mainly controlled by the PIS parameter while the contraction speed was mainly controlled by the DMP parameter.

3.2.4 Delivery of the Fetus.

The simulation results for the delivery of the fetus caused by the contraction of the uterus are shown in Fig. 13. The uterus contracted significantly and relaxed slightly in each contraction cycle, but the overall uterine deformation increased with simulation time. In this version of the model, the fetal head was smaller than the outlet of the uterus—so no stretch to this region of uterine tissue occurred. In addition, the fetal model was simplified as a rigid body, so that no deformation occurred to any portion of the fetus. The contact between the uterus and the fetus was a surface-to-surface contact with no friction defined.

Fig. 13
The delivery of the fetus through the uterus
Fig. 13
The delivery of the fetus through the uterus
Close modal

According to the mechanical process of labor, the movement of the fetus is mostly caused by the pushing force from the active contraction of the uterus with help from a maternal pushing force acting within the mother's abdomen [10]. The current model has not considered the effect of maternal pushing, which would add an effective outward force acting from the fundus. The sole source of the pushing force from the uterus in this simulation was the active contraction forces caused by the contractile fibers inside of the uterine wall. The contraction within the uterus did not happen in all regions simultaneously, but resulted from a contraction wave that initiated from the top region and propagated to the lower part of the uterus, as represented in the model. The fetus started to move at t = 2 s when the inner surface of the uterus started to contact the fetal body. The contraction of the uterus and accumulated uterine deformation in the longitudinal direction through the three contraction cycles pushed the fetus to move downward through the uterus. At t = 180 s, the fetus' head and neck were delivered out of the uterus. The largest element stress and strain was 0.13 MPa and 50.23% in the top of the fundus through the three cycles. The overall displacement of the fetus was around 175 mm at the end of the simulation.

The results of this model simulating the delivery of a fetus from a uterus by the uterine active contraction have demonstrated that bonding the three-dimensional contractile fibers with solid portion of uterine muscle wall can simulate the uterine active contraction successfully. In this case, the uterus contracts actively without the need for any external loading, which aligns with what is observed clinically. The model was able to produce a large enough deformation to push the fetus to move through the uterus.

3.2.5 Validation of Uterine Model.

As is typically true for models related to childbirth, direct validation of this model is difficult as the necessary parameters cannot be safely measured during childbirth. Thus, phenomenological validation and verification were used to confirm that the model is behaving in a manner that matches clinical observations. First, the value of the stress was higher in the fundus region as shown in Fig. 13, which agrees with experimental findings that the intensity of the uterine contraction is stronger in the top region of the uterus [13]. Second, according to clinical data, the time period for the second stage of labor varies between 2 min and 2 h. In this model, the delivery time was 3 min. Finally, the propagation time of the contraction wave was 12 s (first row of Fig. 9), which also agreed with experimental data [13] and previous simulation results [28,29].

4 Discussion

Delivery of the fetus during childbirth is accompanied by not insignificant risk to both the mother and the baby, and the active contraction behaviors of the uterus are vital in regulating such a process. In this paper, we first investigated the active contraction behaviors of a unit of muscle tissue and then developed an FEM model of uterine activation and contraction driven by contractile fibers, and finally simulated the delivery of the fetus from the uterus. The fibers inside of the uterus in the model were three-dimensional and divided into seven regions. This results in a model that is very flexible and can be used to investigate the effects of the strong anisotropy of fiber distribution on the mechanical behavior of the uterus, which is an improvement over a model with two-dimensional, evenly distributed fibers within the uterus [30]. In addition, this model can simulate the uterine contraction wave propagation from the fundus region to the lower parts of the uterus. In the simulation results, the fundus region experienced the highest levels of stress during each contraction, which matched with experimental findings that the contraction intensity is higher in the fundus area than the middle and lower parts of the uterus [13]. The PIS and DMP parameters in the Hill model can significantly influence the mechanical behaviors for both the isolated muscle tissue and the complete uterus. They can be tuned to produce a contraction pattern that agrees with physiological behavior. Unlike many other simulation models in the literature where the movement of the fetus was imposed as a specific trajectory, the delivery of the fetus in this model was caused by the uterine active contractions without any prescribed curves for the fetus to move along.

The existing multiscale electromechanical coupling models that have been used to simulate uterine active contraction are usually mathematically complex, physiologically complicated, and computationally intensive [28,29], and they therefore would need an extremely large amount of computer memory and high computational speeds to predict macrolevel mechanical behavior. Due to these challenges, these models were implemented with highly simplified geometries for the uterus and were not able to simulate the large deformations that are necessary to deliver a fetus. In addition, for these previous models, only the structure of the uterus was modeled, and they did not include a model of the fetus. In contrast, this model can perform large uterine deformation simulation on a complicated geometry and deliver a fetus. In addition, the developed model only took 5 h to complete its calculations, which is a much simpler and less computationally intensive way to realize that effect, and therefore it is an appropriate and efficient way to investigate such an organ's macroscale biomechanical behaviors.

There are of course limitations with this model, as is the case for all physiological simulations. The intrauterine pressure is the hydrostatic pressure of the amniotic fluid between the fetus and the uterus, which contributes to pushing the fetus through the birth canal. In the second stage of labor, the intrauterine peak pressure varies from 50 mmHg to 150 mmHg [40]. Even though the direct contact between the uterine inner wall with the fetus generated the expected pushing effect on the fetus in this model, it is not clear to what degree the increase in intra-uterine pressure contributes to the fetus' delivery. None of the current simulation models, however, have considered the effect of intra-uterine pressure, including this one. In addition, the pelvic structures, including the bony pelvis and pelvic floor muscles, create a strong resistance to the fetus' delivery process. To better explore the pathomechanics of the entire childbirth process, it will be necessary to build a model system that includes both the uterus and the pelvic structures.

In conclusion, an FEM model to simulate uterine cyclic active contraction driven by three dimensional contractile fibers and capable of delivering a fetus was developed in ls-dyna. In this model, the uterus was able to represent the anisotropy of the fiber distribution as well as perform the large deformation, cyclic active contraction, and propagation of the contraction wave needed to achieve fetal delivery. The ability to phenomenologically model the labor and delivery process will increase our understanding of how variations in anatomy or the characteristics of labor might impact injury risk to both the mother and the infant. While not biofidelic on a micro or cellular level, this phenomenological model will be able to demonstrate whether or not the uterus—with variations in its physiological response—can appropriately cause descent of the fetus through the birth canal. This uterine model has been incorporated with pelvic structures so that a model system that more accurately simulates the complete fetus delivery process can be used to investigate the biomechanics of the second stage of labor [41].

Funding Data

  • NSF (Award ID: CBET-2028474; Funder ID: 10.13039/100000001).

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

References

1.
Alrubaii
,
B. J.
,
2005
, “
Lack of Progress in Labour as a Reason for Cesarean
,”
Med. J. Babylon
,
2
(
1
), pp.
589
595
. 10.1016/S0029-7844(99)00575-X
2.
Myers
,
K. M.
, and
Elad
,
D.
,
2017
, “
Biomechanics of the Human Uterus
,”
Wiley Interdiscip. Rev. Syst. Biol. Med.
,
9
(
5
), pp.
1
20
.10.1002/wsbm.1388
3.
Pino
,
J. H.
,
2017
, “
Arrangement of Muscle Fibers in the Myometrium of the Human Uterus: A Mesoscopic Study
,”
MOJ Anat. Physiol.
,
4
(
2
), pp.
131
135
.10.15406/mojap.2017.04.00131
4.
Fujimoto
,
K.
,
Kido
,
A.
,
Okada
,
T.
,
Uchikoshi
,
M.
, and
Togashi
,
K.
,
2013
, “
Diffusion Tensor Imaging (DTI) of the Normal Human Uterus In Vivo at 3 Tesla: Comparison of DTI Parameters in the Different Uterine Layers
,”
J. Magn. Reson. Imaging
,
38
(
6
), pp.
1494
1500
.10.1002/jmri.24114
5.
Fiocchi
,
F.
,
Nocetti
,
L.
,
Siopis
,
E.
,
Currà
,
S.
,
Costi
,
T.
,
Ligabue
,
G.
, and
Torricelli
,
P.
,
2012
, “
In Vivo 3T MR Diffusion Tensor Imaging for Detection of the Fibre Architecture of the Human Uterus: A Feasibility and Quantitative Study
,”
Br. J. Radiol.
,
85
(
1019
), pp.
e1009
e1017
.10.1259/bjr/76693739
6.
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 Discov. Mol., Cellular, Evol. Biol.
,
288A
(
1
), pp.
84
90
.10.1002/ar.a.20274
7.
McLean
,
J. P.
,
Fang
,
S.
,
Gallos
,
G.
,
Myers
,
K. M.
, and
Hendon
,
C. P.
,
2020
, “
Three-Dimensional Collagen Fiber Mapping and Tractography of Human Uterine Tissue Using OCT
,”
Biomed. Opt. Express
,
11
(
10
), p.
5518
.10.1364/BOE.397041
8.
Fang
,
S.
,
McLean
,
J.
,
Shi
,
L.
,
Vink
,
J. S. Y.
,
Hendon
,
C. P.
, and
Myers
,
K. M.
,
2021
, “
Anisotropic Mechanical Properties of the Human Uterus Measured by Spherical Indentation
,”
Ann. Biomed. Eng.
,
49
(
8
), pp.
1923
1942
.10.1007/s10439-021-02769-0
9.
He
,
Y.
,
Ding
,
N.
,
Li
,
Y.
,
Li
,
Z.
,
Xiang
,
Y.
,
Jin
,
Z.
, and
Xue
,
H.
,
2015
, “
3-T Diffusion Tensor Imaging (DTI) of Normal Uterus in Young and Middle-Aged Females During the Menstrual Cycle: Evaluation of the Cyclic Changes of Fractional Anisotropy (FA) and Apparent Diffusion Coefficient (ADC) Values
,”
Br. J. Radiol.
,
88
(
1049
), p.
20150043
.10.1259/bjr.20150043
10.
Grimm
,
M. J.
,
2021
, “
Forces Involved With Labor and Delivery—A Biomechanical Perspective
,”
Ann. Biomed. Eng.
,
49
(
8
), pp.
1819
1835
.10.1007/s10439-020-02718-3
11.
Smith
,
R.
,
Imtiaz
,
M.
,
Banney
,
D.
,
Paul
,
J. W.
, and
Young
,
R. C.
,
2015
, “
Why the Heart is Like an Orchestra and the Uterus is Like a Soccer Crowd
,”
Am. J. Obstet. Gynecol.
,
213
(
2
), pp.
181
185
.10.1016/j.ajog.2015.06.040
12.
Steer
,
P. J.
,
Little
,
D. J.
,
Lewis
,
N. L.
,
Kelly
,
M. C. M. E.
, and
Beard
,
R. W.
,
1975
, “
Uterine Activity in Induced Labour
,”
BJOG
,
82
(
6
), pp.
433
441
.10.1111/j.1471-0528.1975.tb00666.x
13.
Caldeyro-Barcia
,
R.
,
Alvarez
,
H.
, and
Poseiro
,
J. J.
,
1955
, “
Normal and Abnormal Uterine Contractility in Labour
,”
Triangle
,
2
(
41
).
14.
Yao
,
W.
,
Yoshida
,
K.
,
Fernandez
,
M.
,
Vink
,
J.
,
Wapner
,
R. J.
,
Ananth
,
C. V.
,
Oyen
,
M. L.
, and
Myers
,
K. M.
,
2014
, “
Measuring the Compressive Viscoelastic Mechanical Properties of Human Cervical Tissue Using Indentation
,”
J. Mech. Behav. Biomed. Mater.
,
34
, pp.
18
26
.10.1016/j.jmbbm.2014.01.016
15.
Fodera
,
D. M.
,
Russell
,
S. R.
,
Jackson
,
J. L. L.
,
Fang
,
S.
,
Chen
,
X.
,
Vink
,
J.
,
Oyen
,
M. L.
, and
Myers
,
K. M.
,
2023
, “
Material Properties of Nonpregnant and Pregnant Human Uterine Layers
,”
J. Mech. Behav. Biomed. Mater.
,
151
, p.
106348
.10.1016/j.jmbbm.2023.106348
16.
Myers
,
K. M.
,
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
.10.1115/1.3197847
17.
Shi
,
L.
, and
Myers
,
K.
,
2023
, “
A Finite Porous-Viscoelastic Model Capturing Mechanical Behavior of Human Cervix Under Multi-Step Spherical Indentation
,”
J. Mech. Behav. Biomed. Mater.
, 143, p.
105875
.10.1016/j.jmbbm.2023.105875
18.
Shi
,
L.
,
Hu
,
L.
,
Lee
,
N.
,
Fang
,
S.
, and
Myers
,
K.
,
2022
, “
Three-Dimensional Anisotropic Hyperelastic Constitutive Model Describing the Mechanical Response of Human and Mouse Cervix
,”
Acta Biomater.
,
150
, pp.
277
294
.10.1016/j.actbio.2022.07.062
19.
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
.10.1115/1.4036259
20.
Westervelt
,
A. R.
, and
Myers
,
K. M.
,
2017
, “
Computer Modeling Tools to Understand the Causes of Preterm Birth
,”
Semin. Perinatol.
,
41
(
8
), pp.
485
492
.10.1053/j.semperi.2017.08.007
21.
Scott
,
A. K.
,
Louwagie
,
E. M.
,
Myers
,
K. M.
, and
Oyen
,
M. L.
,
2023
, “
Biomechanical Modeling of Cesarean Section Scars and Scar Defects
,” bioRxiv.10.1101/2023.11.03.565565
22.
Koh
,
C. T.
, and
Oyen
,
M. L.
,
2015
, “
Toughening in Electrospun Fibrous Scaffolds
,”
APL Mater
,
3
(
1
), p.
014908
.10.1063/1.4901450
23.
Clark
,
A. R.
,
Yoshida
,
K.
, and
Oyen
,
M. L.
,
2022
, “
Computational Modeling in Pregnancy Biomechanics Research
,”
J. Mech. Behav. Biomed. Mater.
,
128
, p.
105099
.10.1016/j.jmbbm.2022.105099
24.
Parente
,
M. P. L.
,
Natal Jorge
,
R. M.
,
Mascarenhas
,
T.
,
Fernandes
,
A. A.
, and
Martins
,
J. A. C.
,
2009
, “
The Influence of the Material Properties on the Biomechanical Behavior of the Pelvic Floor Muscles During Vaginal Delivery
,”
J. Biomech.
,
42
(
9
), pp.
1301
1306
.10.1016/j.jbiomech.2009.03.011
25.
Xuan
,
R.
,
Yang
,
M.
,
Gao
,
Y.
,
Ren
,
S.
,
Li
,
J.
,
Yang
,
Z.
,
Song
,
Y.
,
Huang
,
X. H.
,
Teo
,
E. C.
,
Zhu
,
J.
, and
Gu
,
Y.
,
2021
, “
A Simulation Analysis of Maternal Pelvic Floor Muscle
,”
Int. J. Environ. Res. Public Health
,
18
(
20
), p.
10821
.10.3390/ijerph182010821
26.
Buttin
,
R.
,
Zara
,
F.
,
Shariat
,
B.
, and
Redarce
,
T.
,
2009
, “
A Biomechanical Model of the Female Reproductive System and the Fetus for the Realization of a Childbirth Virtual Simulator
,”
31st Annual International Conference of the IEEE Engineering in Medicine and Biology Society: Engineering the Future of Biomedicine, EMBC 2009
, Minneapolis, MN, Sept. 2–6, pp.
5263
5266
.10.1109/IEMBS.2009.5334085
27.
Buttin
,
R.
,
Zara
,
F.
,
Shariat
,
B.
,
Redarce
,
T.
,
Buttin
,
R.
,
Zara
,
F.
,
Shariat
,
B.
, et al.,
2013
, “
Biomechanical Simulation of the Fetal Descent Without Imposed Theoretical Trajectory
,”
Comput. Methods. Programs. Biomed.
,
111
(
2
), pp.
389
401
.10.1016/j.cmpb.2013.04.005
28.
Cochran
,
A. L.
, and
Gao
,
Y.
,
2015
, “
A Model and Simulation of Uterine Contractions
,”
Math. Mech. Solids
,
20
(
5
), pp.
540
564
.10.1177/1081286513507940
29.
Yochum
,
M.
,
Laforet
,
J.
, and
Marque
,
C.
,
2016
, “
An Electro-Mechanical Multiscale Model of Uterine Pregnancy Contraction
,”
Comput. Biol. Med.
,
77
, pp.
182
194
.10.1016/j.compbiomed.2016.08.001
30.
Vila Pouca
,
M. C. P.
,
Ferreira
,
J. P. S.
,
Oliveira
,
D. A.
,
Parente
,
M. P. L.
,
Mascarenhas
,
M. T.
, and
Natal Jorge
,
R. M.
,
2019
, “
Simulation of the Uterine Contractions and Foetus Expulsion Using a Chemo-Mechanical Constitutive Model
,”
Biomech. Model. Mechanobiol.
,
18
(
3
), pp.
829
843
.10.1007/s10237-019-01117-5
31.
Hill
,
A. V.
,
1938
, “
The Heat of Shortening and the Dynamic Constants of Muscle
,”
Proc. R. Soc. London. Ser. B
,
126
(
843
), pp.
136
195
.10.1098/rspb.1938.0050
32.
Bates
,
J. H. T.
, and
Lauzon
,
A.-M.
,
2007
, “
Parenchymal Tethering, Airway Wall Stiffness, and the Dynamics of Bronchoconstriction
,”
J. Appl. Physiol.
,
102
(
5
), pp.
1912
1920
.10.1152/japplphysiol.00980.2006
33.
Murphy
,
R. A.
,
1988
, “
Muscle Cells of Hollow Organs
,”
Physiology
,
3
(
3
), pp.
124
128
.10.1152/physiologyonline.1988.3.3.124
34.
Östh
,
J.
,
2014
,
Muscle Responses of Car Occupants: Numerical Modeling and Volunteer Experiments Under Pre-Crash Braking Conditions
,
Chalmers Tekniska Hogskola
,
Sweden
.
35.
Bamberg
,
C.
,
Rademacher
,
G.
,
Güttler
,
F.
,
Teichgräber
,
U.
,
Cremer
,
M.
,
Bührer
,
C.
,
Spies
,
C.
, et al.,
2012
, “
Human Birth Observed in Real-Time Open Magnetic Resonance Imaging
,”
Am. J. Obstet. Gynecol.
,
206
(
6
), pp.
505-e1
505.e6
.10.1016/j.ajog.2012.01.011
36.
Zhang
,
W.
, and
Chen
,
J.
,
2020
, “
Diffusion Tensor Imaging (DTI) of the Cesarean-Scarred Uterus In Vivo at 3T: Comparison Study of DTI Parameters Between Nonpregnant and Pregnant Cases
,”
J. Magn. Reson. Imaging
,
51
(
1
), pp.
124
130
.10.1002/jmri.26868
37.
Dubrauszky
,
V.
,
Schwalm
,
H.
, and
Fleischer
,
M.
,
1971
, “
The Fibre System of Connective Tissue in the Childbearing Age, Menopause, and Pregnancy
,”
Arch. Gynakol.
,
210
(
3
), pp.
276
292
.10.1007/BF00667740
38.
Östh
,
J.
,
Brolin
,
K.
, and
Bråse
,
D.
,
2015
, “
A Human Body Model With Active Muscles for Simulation of Pretensioned Restraints in Autonomous Braking Interventions
,”
Traffic Inj. Prev.
,
16
(
3
), pp.
304
313
.10.1080/15389588.2014.931949
39.
Abalos
,
E.
,
Oladapo
,
O. T.
,
Chamillard
,
M.
,
Díaz
,
V.
,
Pasquale
,
J.
,
Bonet
,
M.
,
Souza
,
J. P.
, and
Gülmezoglu
,
A. M.
,
2018
, “
Duration of Spontaneous Labour in ‘Low-Risk’ Women With ‘Normal’ Perinatal Outcomes: A Systematic Review
,”
Eur. J. Obstet. Gynecol. Reprod. Biol.
,
223
, pp.
123
132
.10.1016/j.ejogrb.2018.02.026
40.
Lindgren
,
L.
,
1960
, “
The Causes of Foetal Head Moulding in Labour
,”
Acta Obstet. Gynecol. Scand.
,
39
(
1
), pp.
46
62
.10.3109/00016346009157836
41.
Tao
,
R.
, and
Grimm
,
M. J.
,
2024
, “
Simulation of the Childbirth Process in LS-DYNA
,”
ASME J. Biomech. Eng.
,
146
(
6
), p.
061002
.10.1115/1.4064594