Biomech Model Mechanobiol DOI 10.1007/s10237-015-0703-z

ORIGINAL PAPER

Simulating uterine contraction by using an electro-chemo-mechanical model Babak Sharifimajd1 · Carl-Johan Thore1 · Jonas Stålhand1

Received: 8 January 2015 / Accepted: 2 July 2015 © Springer-Verlag Berlin Heidelberg 2015

Abstract Contractions of uterine smooth muscle cells consist of a chain of physiological processes. These contractions provide the required force to expel the fetus from the uterus. The inclusion of these physiological processes is, therefore, imperative when studying uterine contractions. In this study, an electro-chemo-mechanical model to replicate the excitation, activation, and contraction of uterine smooth muscle cells is developed. The presented modeling strategy enables efficient integration of knowledge about physiological processes at the cellular level to the organ level. The model is implemented in a three-dimensional finite element setting to simulate uterus contraction during labor in response to electrical discharges generated by pacemaker cells and propagated within the myometrium via gap junctions. Important clinical factors, such as uterine electrical activity and intrauterine pressure, are predicted using this simulation. The predictions are in agreement with clinically measured data reported in the literature. A parameter study is also carried out to investigate the impact of physiologically related parameters on the uterine contractility. Keywords Excitation–contraction model of uterine smooth muscle · Uterus contraction · Intrauterine pressure · Uterine electrical activity

B 1

Babak Sharifimajd [email protected] Division of Solid Mechanics, Department of Management and Engineering, The Institute of Technology, Linköping University, 581 83 Linköping, Sweden

1 Introduction Complications during labor are mainly due to uterine contraction difficulties (Gifford et al. 2000; Maul et al. 2003; Zhu et al. 2006). Preterm labor, for example, is caused by early onset of uterine contractions and causes 35 % of all neonatal deaths worldwide (Blencowe et al. 2012), as well as significant levels of long-term morbidity, including neurodevelopmental deficits and increased chronic disease risk in adulthood (Mwaniki et al. 2012). Further, prolonged labors which occur in 10 % of all pregnancies and lead to 20 % of cesarean sections (Bugg et al. 2006) are also due to dysfunctional uterine contraction (Johnson and Clayton 1955). An improved understanding of the physiological complexities of myometrial activity will, therefore, assist in the task of developing better-targeted therapies for these problematic labors. Ethical considerations, however, usually prevent direct studies on the human uterus. This leaves us no choice but to develop models to simulate uterine contractions and study the effects of different parameters involved in this process on the uterine contractility. Uterine smooth muscle cells (USMCs), the main constituent of the myometrium, generate uterine contractions through a sequence of coupled physiological processes, referred to as the excitation–contraction (EC) process (Sandow 1952). Through this process, the cell membrane gets excited by external (electrical) stimuli, and, in response, the transmembrane (electrical) potential undergoes a rapid change, known as an action potential (AP) that opens membrane channels to let ions enter or leave the cell (Keener and Sneyd 2008). Among these ions, calcium ions (Ca2+ ) are important since they trigger myosin phosphorylation process which leads to cross-bridge formation and cycling, and eventually cell contraction (Hai and Murphy 1988). The coupling

123

B. Sharifimajd et al.

between these processes requires an integrated model which can replicate the muscle response in a physiologically consistent manner. Surprisingly, there are no studies that investigate uterine contraction based on the EC coupling process of USMCs at an organ level. Most of the published studies on uterine activity concern only the electrophysiology of the USMC membrane (Buhimschi et al. 1997; Maul et al. 2003; Burdyga et al. 2007; Bursztyn et al. 2007; Tong et al. 2011; La Rosa et al. 2012). The only study which has simulated uterine contraction is due to Cochran and Gao (2013). They developed an electromechanical model to simulate intrauterine pressure (IUP), a clinical gold standard to assess the uterine contraction during labor (Csapo and Pinto-Dantas 1965; Maul et al. 2003; Buhimschi et al. 2004; Bastos et al. 2010). The model does not, however, include important parts of the process leading to uterine contraction, namely the dynamics of the intracellular calcium ion concentration ([Ca2+ ]i ), the myosin phosphorylation process, and the filament sliding, known to be the origin of muscle contractions (Huxley 1974; McMahon 1984; Murtada et al. 2010; Sharifimajd and Stålhand 2013). Instead it relates the active stress directly to the AP through a single, phenomenological differential equation. This causes the active stress behavior to follow the AP shape closely, leading to an IUP which does not resemble clinical measurements. Based on the framework established in Sharifimajd and Stålhand (2014), an electro-chemo-mechanical (ECM) model is developed to replicate the USMC excitation– contraction in this work. Uterine contraction during labor is simulated by implementing this model in the finite element (FE) solver Comsol Multiphysics . Important clinical factors such as IUP and uterine electrical activity are predicted through the simulation. The ECM model proposed herein makes it possible to account for physiological phenomena controlling uterine contractility and study their effects on the IUP.

2 The electro-chemo-mechanical model The ECM model consists of three sub-models: – A membrane model which includes a FitzHugh–Nagumo type model (FHN) (FitzHugh 1961; Nagumo et al. 1962; Aliev and Panfilov 1996) to simulate AP coupled to a model that determines the [Ca2+ ]i dynamics. – A phosphorylation model to determine the fraction of bound myosin (cross-bridges). – A mechanical model to compute passive and active responses of the uterus. In what follows, these sub-models will be presented in detail.

123

2.1 Membrane model: electrical activity The smooth muscle cell membrane is a very complex structure with many channels and gates for transporting ions to and from the intracellular space. The membrane can be triggered by electrical signals (impulses from neurons or pacemaker cells), chemical agents (hormones), or mechanical means (deformation). The overall effect is an increase in the [Ca2+ ]i which leads to muscle contraction. 2.1.1 Membrane excitability and AP generation The process which leads to AP generation can be replicated through a Hodgkin–Huxley-like model, see, e.g., (Hodgkin and Huxley 1952; Burdyga et al. 2007; Tong et al. 2011). In Sharifimajd and Stålhand (2014), it is shown that the ionic currents can be determined constitutively by defining an electrochemical potential. The total current defined in the Hodgkin–Huxley model then follows from conservation of charge for all ions entering and leaving the cell. The Hodgkin–Huxley model has been vastly used to explain the excitability of smooth muscle cell membranes (Victorri et al. 1985; Tong et al. 2011; Yang et al. 2003; Toth et al. 2013; Korogod and Kochenov 2013; Sharifimajd and Stålhand 2014). The model is complex and computationally expensive, however. Many studies, including the present, therefore, use a reaction–diffusion equation of FHN type to mimic AP generation and diffusion (Benson et al. 2006; Göktepe and Kuhl 2010; Ambrosi et al. 2011; Nobile et al. 2012; La Rosa et al. 2012; Cochran and Gao 2013). We start by introducing U0 and U as the regions in space occupied by the uterus at time t = 0 and t > 0, respectively. The former represents the reference configuration, while the latter represents the current configuration. Each point X in U0 is mapped to a point x = x(X, t) in U. Following Aliev and Panfilov (1996), the reaction–diffusion equation is taken to be ∂t u + ∇ · (u x˙ ) = ∇ · (d∇u) in U, + cu(u − a)(1 − u) − r u + Ist ,   in U, ∂t r + ∇ · (r x˙ ) = b cu(1 + a − u) − r ,

(1a) (1b)

where ∂t and a superscribed dot indicates, respectively, a spatial and a material time derivative. The scalar field u(x, t) is a voltage-like variable, and r (x, t) is the recovery variable (FitzHugh 1961; Aliev and Panfilov 1996), both defined per unit volume. The second-order tensor d is a spatial diffusion tensor; Ist is the external, current-like stimulus; and a, b, c are positive constants. The symbols ∇ and ∇· denote gradient and divergence operators in spatial coordinates, respectively. The divergence terms in the left-hand side of Eqs. (1a) and (1b) indicate convective transport due to the displacement

Simulating uterine contraction by using an electro-chemo-mechanical model

of the material itself. Since x˙ is small during smooth muscle contraction, these terms become very small and can be neglected (Ambrosi et al. 2011). For a material description of Eqs. (1a) and (1b), see Remark 1. Equations (1a) and (1b) are complemented with the boundary condition d∇x u · n = 0,

on ∂U,

(2)

where n is the spatial normal vector to the boundary ∂U. The values of u will not fall into the physiological range of an AP for smooth muscle cells (Keener and Sneyd 2008). A simple step function is therefore defined to compute the transmembrane potential Vm = Vm (X, t) from u(x(X, t), t) as  res Vm if u < u, ¯ Vm = (3) dep Vm otherwise, dep

where Vmres is the resting potential, Vm is the depolarization potential, and u¯ ∈ [u min , u max ] with u min and u max being the minimum and maximum values of u, respectively. 2.1.2 Intracellular calcium ion dynamics The FHN-type model will not provide information regarding the dynamics of [Ca2+ ]i . The simple model developed by Bursztyn et al. (2007) has been, therefore, chosen to determine the [Ca2+ ]i from the APs. Even though the model is simple, it recovers the main features of [Ca2+ ]i raise in USMCs after membrane depolarization by assuming: Ca2+

entry following membrane depolarization occurs almost entirely via L-type voltage-gated Ca2+ channels (VOCCs); and – Ca2+ extraction out of the cell occurs by Ca2+ -pumps and Na+ /Ca2+ -exchangers in the membrane. –

The calcium influx via L-type VOCCs can be calculated by using JVOCC

gCa = (E Ca − Vm ) , z Ca F Vcell

(4)

where gCa is the conductivity of L-type VOCCs to be specified later, Vcell is the cell volume, and E Ca is the Nernst potential for Ca2+ defined by E Ca =

  [Ca2+ ]e RT , ln z Ca F [Ca2+ ]i

(5)

where R is the ideal gas constant, T is the absolute temperature, F is the Faraday constant, z Ca is the valence of Ca2+ , and [Ca2+ ]e is the extracellular calcium concentration, assumed

constant herein. The Ca2+ conductivity gCa is described by a Boltzmann-type activation curve as follows (Bursztyn et al. 2007; Parkington et al. 1999): gCa =

1 + exp



gm,Ca  , VCa,1/2 − Vm /K Ca,1/2

(6)

where gm,Ca is the maximum Ca2+ conductance, VCa,1/2 is the half-activation potential of L-type VOCCs, and K Ca,1/2 is the slope of the activation function at VCa,1/2 (Hodgkin and Huxley 1952). The efflux of Ca2+ via the Ca2+ -pumps can be described by the Hill function n

Jpump = Vpmax

[Ca2+ ]i , [Kph ]n + [Ca2+ ]i n

(7)

where n is the Hill coefficient, K ph is the Ca2+ concentration in half-activation of the pump, and Vpmax is the maximal rate of Ca2+ extraction by the pump (Bursztyn et al. 2007; Keener and Sneyd 2008). The last component involved in [Ca2+ ]i dynamics is the Na+ /Ca2+ -exchanger which can extract calcium from the cell by substituting Ca2+ for extracellular Na+ (Parthimos et al. 1999; Keener and Sneyd 2008). The flux of Ca2+ through all Na+ /Ca2+ -exchangers is determined by (Bursztyn et al. 2007) JNa/Ca = gNa/Ca

  [Ca2+ ]i Vm − E Na/Ca , 2+ [Ca ]i + K Na/Ca

(8)

where gNa/Ca is the conductance of all Na+ /Ca2+ -exchangers within the cell membrane, K Na/Ca is the intracellular Ca2+ concentration in half-activation of the Na+ /Ca2+ exchangers, and E Na/Ca is the reversal potential of the exchangers given by E Na/Ca = 3E Na − 2E Ca ,

(9)

where E Na is the Nernst potential for sodium defined by E Na

  [Na+ ]e RT , ln = F [Na+ ]i

(10)

where [Na+ ]e and [Na+ ]i are constant extracellular and intracellular Na+ concentrations, respectively. The rate-of-change of [Ca2+ ]i is now computed through ˙ [Ca2+ ]i = JVOCC − Jpump + JNa/Ca

in U0 .

(11)

An ionic influx from extracellular into intracellular space is considered to be negative. In addition, the transient effects of Ca2+ buffering by the superficial SR are not taken into account in this model (Bursztyn et al. 2007).

123

B. Sharifimajd et al. Fig. 1 A schematic illustration of myosin phosphorylation kinetics

2.3 Uterine smooth muscle contraction model: mechanical activity

2.2 Myosin phosphorylation model: chemical activity It is known that Ca2+ inside smooth muscle cells can bind to calmodulin (CaM) and form Ca2+ -CaM complex which activates myosin light-chain kinase (MLCK). MLCK phosphorylates the regulatory light chains of the myosin heads. Phosphorylated myosin heads are then able to bind actin filaments and form cross-bridges. This indicates that the activation level of smooth muscle, unlike other muscles, depends on the level of phosphorylation of myosin heads and not the applied stimulus directly. The kinetics of myosin phosphorylation, herein, is taken to be described by the model in Hai and Murphy (1988) which assumes that myosin can be in one of four states (Fig. 1): unbound unphosphorylated (α A ), unbound phosphorylated (α B ), bound phosphorylated (αC ), and bound dephosphorylated (α D ). Since the states are fractions, they must satisfy the constraint α A +α B +αC +α D = 1. The rates-of-change between the four myosin states is given by ⎤ ⎡ ⎤ −k1 α˙ A k2 0 k7 ⎢ α˙ B ⎥ ⎢ k1 −(k2 + k3 ) ⎥ k4 0 ⎢ ⎥=⎢ ⎥ ⎣ α˙ C ⎦ ⎣ 0 ⎦ −(k4 + k5 ) k6 k3 α˙ D 0 0 k5 −(k6 + k7 ) ⎡ ⎤ αA ⎢α B ⎥ ⎥ ×⎢ in U0 , (12) ⎣ αC ⎦ αD ⎡

where k1 to k7 are reaction rates. Following Hai and Murphy (1988) and Bursztyn et al. (2007), the phosphorylation rate k1 (X, t) is taken to be dependent on [Ca2+ ]i through

k1 =

[Ca2+ ]i [Ca2+ ]i

n k1

n k1

+ [Ca2+ (1/2MLCK) ]i

n k1

,

(13)

where [Ca2+ (1/2MLCK) ]i is the required intracellular calcium ion concentration for half-activation of MLCK by the calcium–calmodulin complex, and n k1 is the Hill coefficient of activation. Furthermore, it is assumed that k6 = k1 , k2 = k5 , and k3 = 4k4 (Hai and Murphy 1988; Bursztyn et al. 2007).

123

A continuum model for smooth muscle contraction was derived from basic principles of thermodynamics and mechanics in Sharifimajd and Stålhand (2014). This model is extended herein to determine the mechanics of the uterus including different uterine smooth muscle (USM) fiber families and other passive structures. The interested reader is referred to the original paper for a complete derivation of the model. 2.3.1 Filament sliding kinematics Smooth muscle contraction results from a relative sliding between actin and myosin filaments. Following Sharifimajd and Stålhand (2014), a contractile element (CE) is defined to mimic the dynamics of the contractile machinery within USMC and represented by the mechanical model shown in Fig. 2. The overall effect of cycling cross-bridges is modeled by a friction clutch which displaces actin by applying a friction force FC trough its rotation. There is a brakelike mechanism in series with the clutch which represents the overall effect of all non-cycling, or slow-cycling, crossbridges, also known as latch bridges (Fig. 1). The clutch and the brake are attached to an elastic element which represents the elasticity of all cross-bridges within the CE. The CE is assumed to deform from its initial length l0 to the final length l in two steps: first, by an active displacement u¯ a in which only the clutch and brake are displaced relative to actin while the elastic element is unstretched; and, second, by an elastic displacement u¯ e of the elastic element while the clutch and brake remain fixed relative to actin, see Fig. 2. The uniaxial stretch λ¯ = l/l0 of the CE can then be described by λ¯ = 1 + ε¯ a + ε¯ e ,

(14)

where ε¯ a = u¯ a /l0 and ε¯ e = u¯ e /l0 denote the active and elastic strain, respectively. 2.3.2 Uterine smooth muscle tissue contraction kinematics The myometrium is assumed to comprise two USM families of fasciculi which are parallel to each other and to the surface of the uterus (Young and Hession 1999). These USM families are taken to be aligned with unit vectors Mi (i = 1, 2) defined in the reference configuration. A parameter ζ1 ∈ [0, 1] is introduced to represent the fraction of all fasciculi that belong to the first family. The deformation gradient is defined by F = ∂x/∂X and decomposed into a volumetric part J 1/3 I, and an isochoric (volume preserving) part F¯ = J −1/3 F, with I being the second-order identity tensor and J = det F > 0 the volume

Simulating uterine contraction by using an electro-chemo-mechanical model

Fig. 2 An illustration of the multi-scale uterine smooth muscle histology. From left to right: pregnant uterus [used with the permissions from Londino (2014)], uterine smooth muscle fascicle, USMCs, smooth muscle cell cytoplasm, contractile unit, and the mechanical model for the contractile unit

ratio. It is assumed that the USMCs within each family are structurally homogeneous. Hence, the CE stretch λ¯ i within the USM family i, in the direction Mi , can be described ¯ = F¯ T F¯ through the isochoric right Cauchy–Green tensor C as ¯ : (Mi ⊗ Mi ), λ¯ i2 = C

To describe the mechanical behavior of the uterus, we postulate that there exists a strain-energy function Ψ which can be decoupled into an isochoric contribution Ψ¯ and a volumetric contribution U according to Ψ = Ψ¯ + U.

(17)

(15)

where: denotes double contraction and ⊗ the dyadic product (Holzapfel 2000). 2.3.3 Mechanical balance laws and constitutive relations In the absence of body forces and by neglecting the inertia, the uterine contraction mechanics is governed by the balance of linear momentum together with boundary and initial conditions: ∇x · σ = 0, σ n = t¯,

in U,

(16a)

on ∂Ut ,

(16b)

¯ u = u,

on ∂Uu ,

(16c)

u(x, t)|t=0 = 0,

in U,

(16d)

where σ (x, t) is the Cauchy stress tensor, t¯ is the prescribed surface traction on the surface ∂Ut , u(x, t) = x − X(x, t) is the spatial displacement field, and u¯ is the prescribed dis placement of the surface ∂Uu . Note that ∂Uu ∂Ut = ∂U and ∂Uu ∂Ut = ∅. The uterus has three layers: endometrium, myometrium, and perimetrium (Snell 2011). In this study, however, the uterus is assumed to comprise one homogeneous layer which reflects the mechanical properties of all three layers, without making explicit the contribution of the endo- and perimetrium to the passive mechanical properties of the wall.

The volumetric part of the strain energy is taken to be (Simo and Taylor 1991) U=

 κvol  2 J − 1 − 2 ln(J ) , 4

(18)

where κvol > 0 is a constant. Inspired by the mechanics of fiber-reinforced composites (Spencer 1984) and its applications to arterial wall mechanics (Holzapfel et al. 2000; Böl et al. 2012; Murtada and Holzapfel 2014), and USMC mechanics (Cochran and Gao 2013) together with the mechanical model presented in Sharifimajd and Stålhand (2014), the isochoric strain-energy function is taken to be   p Ψ¯ = Ψ¯ m ( I¯1 ) + ζ1 Ψ¯ p ( I¯4 ) + N1 (¯εa1 )Ψ¯ a (¯εe1 , αxb )   + (1 − ζ1 ) Ψ¯ p ( I¯6 ) + N2 (¯εa2 )Ψ¯ a (¯εe2 , αxb ) , (19) where Ψ¯ m represents the energy stored in the isotropic passive matrix, Ψ¯ p represents the energy stored in passive USM fasciculi and collagen fibers, and Ψ¯ a (active strain energy) represents the elastic energy stored in cross-bridges. The functions Ni : R → [0, 1] account for the overlap between actin and myosin filaments (Stålhand et al. 2011; Sharifimajd and Stålhand 2013, 2014) where the unity value corresponds to the maximum overlap between filaments within the CE. The fraction of bound myosin is given by αxb = αC + α D , p

123

B. Sharifimajd et al.

and ε¯ ai and ε¯ ei are, respectively, the active and elastic strains in the ith USM family. Finally, ¯ : I, ¯ : (M1 ⊗ M1 ), I¯4 = C I¯1 = C ¯ : (M2 ⊗ M2 ), I¯6 = C

(20)

are the first, fourth, and sixth invariant of the isochoric Cauchy–Green tensor, respectively (Holzapfel et al. 2000). Following Holzapfel et al. (2000) and Cochran and Gao (2013), it is assumed that the behavior of the passive matrix can be described by p Ψ¯ m =

 μp  I¯1 − 3 , 2

(21)

where μ p > 0 is a constant. The strain energy Ψ¯ p is determined by (Holzapfel et al. 2000) Ψ¯ p =



    (c1 /2c2 ) exp c2 ( I¯j − 1)2 − 1 0

if I¯j ≥ 1 otherwise,

where ε˙¯ ai = ε˙¯ ai (X, t) is the active strain-rate, η > 0 is a constant, and FC,i = FC,i (X, t) is the friction force generated by the clutch in the ith fiber family (Sharifimajd and Stålhand 2014). This force is related to the cycling cross-bridges and can be described through  FC,i =

Ni κC αC (ε˙¯ ai + vC ) 0

if ε˙¯ ai > −vC , otherwise,

(27)

where κC > 0 is a constant, and vC is the peripheral clutch velocity, see Fig. 2. Since muscle contraction is due to filament sliding, the active force is considered to be dependent on the amount of overlap between actin and myosin. The overlap between filaments is described by (Sharifimajd and Stålhand 2013, 2014) 2 ⎞  opt − ε¯ a − ε¯ ai ⎟ ⎜ Ni = exp ⎝ ⎠, 2γ 2 ⎛

(28)

(22) opt

where j ∈ {4, 6}, and c1 and c2 are positive constants. The elastic energy stored in the cross-bridges for the ith USM family is taken to be (Sharifimajd and Stålhand 2013, 2014) Ψ¯ a =

1 2 αxb κxb ε¯ ei , 2

(23)

where κxb > 0 is a constant. The Cauchy stress tensor is decoupled according to (Holzapfel et al. 2000) σ = σ¯ + σvol ,

(24)

where σ¯ is the isochoric Cauchy stress and σvol is the volumetric Cauchy stress, defined by  ∂U I, σvol = ∂J   ¯ −1 ¯ ∂ Ψ ¯ T σ¯ = P : 2J F F , ¯ ∂C

where ε¯ a denotes the active strain at which the muscle generates maximum force, and γ is a constant. This form of Ni reflects that length dependence of the active force in the smooth muscle is known to be bell-shaped (Herlihy and Murphy 1973; Uvelius 1976; Arner 1982). Finally, to define the diffusion tensor introduced in Eq. (1a), it is assumed that an AP can diffuse anisotropically along fasciculi and isotropically through the matrix. A parameter ξiso ∈ [0, 1] is introduced to represent the fraction of isotropic diffusion. Since the elongation of USMC fibers will not affect AP propagation (Böl et al. 2012), the anisotropic part of the diffusion is taken to depend only on the rotation of the USM fibers through the rotated structural tensors miR defined by miR = R(Mi ⊗ Mi )RT ,

(29)



(25a) (25b)

where P = I − (I ⊗ I)/3 denotes the fourth-order spatial projection tensor including the fourth-order identity tensor I (Holzapfel 2000). An evolution law to describe the dynamics of the friction clutch is derived in Sharifimajd and Stålhand (2014). For USM family i, this law reads ηε˙¯ ai

∂ Ψ¯ ∂ Ψ¯ = −FC,i − + , ∂ ε¯ ai ∂ ε¯ ei

123

(26)

where R is the rotation tensor obtained from the polar decomposition of the deformation gradient tensor F. The diffusion tensor is defined by   d = ξiso I + (1 − ξiso ) ζ1 m1R + (1 − ζ1 )m2R ,

(30)

where > 0 is the diffusion coefficient, assumed to be the same for both USMCs and the matrix (Cochran and Gao 2013). Remark 1 The reaction–diffusion system in Eqs. (1a) and (1b) can be written in the reference configuration as (Ambrosi et al. 2011)

Simulating uterine contraction by using an electro-chemo-mechanical model

  (J˙u) = Div J D Grad(u)   + J cu(u − a)(1 − u) − r + Ist   (J˙ r ) = J b cu(1 + a − u) − r in U0 ,

in U0 , (31a) (31b)

where Div and Grad denote the divergence and gradient operators in material coordinates, respectively, and the second-order tensor D = F−1 dF−T is referred to as the material (referential) diffusion tensor.

3 Finite element implementation The ECM model presented in Sect. 2 is used to simulate uterine contractions and predict the IUP in response to electrical stimuli which depolarize the USMC membranes. The simulation is performed using the FE solver Comsol Multiphysics . The model equations were solved by using an implicit time-stepping scheme with variable step-size. A ‘General Form PDE’ module for solving the variables u and r , and three ‘Domain ODEs and DAEs’ modules for solving [Ca2+ ]i , α A , α B , αC , α D , and ε¯ ai were coupled to a ‘Structural Mechanics’ module in which the domain was considered to consist of a hyperelastic material, defined by the strain-energy function (17). Table 1 shows the initial values for variables solved in Comsol Multiphysics . Nodal displacements in the FE simulation are initialized to zero; cf. Eq. (16d).

Fig. 3 Geometry of the upper segment of a pregnant uterus

which is between 0.6–1.0 cm at 37 weeks of gestation (Usher and McLean 1969; Buhimschi et al. 2003). As mentioned in Sect. 2.3.2, the uterus is considered to consist of two USM families. Fasciculi in the first family are assumed to be aligned longitudinally with respect to the uterus, while those in the other family are taken to be aligned circumferentially (Fig. 4). This is based on the findings that in non-pregnant human uteri, there are collagen fibers aligned both circumferentially and longitudinally, and collagen fibers are believed to align in the same direction as USM fasciculi (Sobotta 1891; Kawarabayashi and Marshall 1981; Weiss et al. 2006). The material unit vector Mi in a spherical coordinate system (Fig. 4) is given by ⎤ ⎡ cos θi sin α cos ϕ Mi = ⎣ cos θi sin α sin ϕ ⎦ , sin θi sin α

(32)

3.1 Uterus geometry In this study, the upper segment of uterus is only considered whose geometry is taken to be presented by an ellipsoid (Fig. 3). The dimensions of the ellipsoid are taken to be close to the human uterus at 36 weeks of gestation (Danforth 1971; Padubidri and Anand 2006). The thickness of the wall is assumed constant and is set to 0.9 cm. This value is in the range of the experimentally measured myometrial thickness

where θi is the in-wall dispersion of fasciculi, α is the azimuthal angle, and ϕ is the polar angle of the uterus surface; see Fig. 4.

3.2 Boundary conditions The upper segment of the uterus is subject to the following boundary conditions (Fig. 3):

Table 1 The initial values of variables Parameter

Initial value

Equations

u

0

(1a) (1b)

r

−0.1

(1a) (1b)

[Ca2+ ]i

180 nM

(11)

αA

1

(12)

αB

0

(12)

αC

0

(12)

αD

0

(12)

ε¯ ai (i = 1, 2)

0

(26)

– The cervical end is fixed in all directions. – It is assumed that the uterus cavity is filled with an incompressible liquid. – A constant pressure is applied to the outer wall to compensate for the baseline intrauterine pressure (uterine tone) (Buhimschi et al. 2004). – The outer wall of the domain is free of any flux of u; cf. Eq. (2). – The geometry is considered to be symmetric about the Y Z plane.

123

B. Sharifimajd et al. Fig. 4 The USM families aligned longitudinally and circumferentially. The alignment of USM fibers is defined by using the illustrated coordinate system in which the red line indicates the USM fasciculi, with et being the tangential unit vector on the uterus surface

4 Intrauterine pressure simulation Intrauterine pressure is considered as the output of the FE simulation, herein. The IUP has been previously computed as a function of uterus cavity volume through a hydraulic model in Rumpel et al. (2005) and Cochran and Gao (2013). In contrast to these studies, we simulate IUP by taking the uterus cavity to be filled with an incompressible liquid which will be pressurized during contraction. In Comsol Multiphysics , this is modeled by applying a variable pressure to the inner wall of the uterus to keep its volume unchanged. Inspired by the functionality of pacemaker cells in myometrium (La Rosa et al. 2012; Buhimschi et al. 2003), the FE simulation starts by applying Ist to the upper-right side of the domain such that u is generated at that region, causing a local depolarization while the rest of the USMCs are at their resting potential. As Vm propagates through the domain, the USMC membranes are depolarized and [Ca2+ ]i increases. Consequently, the fraction of phosphorylated myosin that can undergo cross-bridge cycling is raised. This will cause the USMCs and uterus to contract, and an IUP to be generated. 4.1 Baseline simulation A plateau-type AP is modeled using Eqs. (1a), (1b), and (3), with the resting and depolarization potentials set to, respectively, −80 and 0 mV (Bursztyn et al. 2007; Tong et al. 2011). The duration of the AP is controlled by Ist , see Fig. 5. The AP diffuses through the uterus as shown in Fig. 6 with a diffusion (propagation) speed of about 4 cm/s which agrees with clinical observations (Rabotti et al. 2010; Lammers 2013).

123

Fig. 5 The AP signal which is determined by the variable u. The stimulus current is applied for 90 s

The simulated IUP is compared to the clinical data in Buhimschi et al. (1997) in Fig. 7a. The pressure applied to the outer boundary of the uterus is set to 28 mmHg to account for the uterine tone (Buhimschi et al. 1997). The IUP is naturally correlated with the active stress generated by USMCs. To obtain the active stress herein, first the total isochoric stress is computed through back-substitution of Eqs. (21)–(23) into (19) and then into (25b). Second, the p passive isochoric stresses, related to Ψ¯ m and Ψ¯ p are subtracted from the total isochoric stress, and the result is the active isochoric stress (Sharifimajd and Stålhand 2014). Figure 7b illustrates the uniaxial active stresses generated by each USM family. By comparing with Fig. 7a, it is seen that the shape of IUP waveform is similar to the active stress waveform. The dynamics of [Ca2+ ]i and myosin states are evaluated at a point in the lower-left side of the domain, as shown in Fig. 8a. The deformation of the uterus is quantified by calculating the total displacement within the domain, as shown

Simulating uterine contraction by using an electro-chemo-mechanical model

Fig. 6 The propagation of u through the uterus as time progress

Fig. 7 a The simulated IUP (solid line) compared to the clinical measurement from Buhimschi et al. (1997) (circles). b The active stress generated within two USM fasciculi

Fig. 8 a Bottom the [Ca2+ ]i dynamics. Top myosin states dynamics. b Total deformation in cm as time progresses

in Fig. 8b. The uterine contraction starts after 10 s, while the total displacement reaches its maximum at t = 100 s. The parameters used in this simulation are either taken from literature or tuned such that an acceptable fit between the simulated IUP, and the clinical data are achieved (Table 2).

4.2 Parameter study The presented ECM model accounts for coupled excitation, activation, and contraction processes of USMCs. Parameters controlling each of these processes can affect the IUP

123

B. Sharifimajd et al. Table 2 The parameters used in the simulation of uterus contraction

Parameter

Value

Equations

Parameter

Value

a

0.003

(1a)

c1

6 kPa

(22)

b

0.02

(1b)

c2

1

(22)

c

5

(1a)

μ†p

60 kPa

(21)

5.729E−5 mM/(mV s)

(8)

κvol

1 MPa

(18)

(30)

κxb

6 MPa

(23)

(19) (30)

ε¯ a

opt

0.4950

(28)

∗ gNa/Ca

†

ζ1†

0.6

mm2 /ms

0.4

Equations

† ξiso

1/11

(30)

γ

1

(28)

R

8314 mJ/(mol K)

(5) (10)

η

18 MPa s

(26)

T

310 K

(5) (10)

vC

0.11 s−1

(27)

F

96,485 C/mol

(5) (10)

κC

3.5 MPa s

(27)

z Ca

2

(5) (4)

k2∗

1.2387 s−1

(12)

0.1419 s−1

(12)

0.0378 s−1

(12)

256.98 nM

(13)

8.7613

(13)

2 mM

(5)

2.98 mM

(10)

140 mM

(10)

∗ Vcell ∗ VCa,1/2 ∗ K Ca,1/2 ∗ K Na/Ca ∗ K ph ∗ gm,Ca ∗ Vpmax ∗ n

21E−15

m3

(4)

−27 mV

(6)

11 mV

(6)

7 µM

(7)

0.6 µM

(7)

0.046842 nS

(6)

0.5145E−3 mM/(s mV)

(7)

1.9015

(7)

k3∗ k7∗

∗ [Ca2+ (1/2MLCK) ]i n ∗k1 [Ca2+ ]∗ext [Na+ ]∗i [Na+ ]∗ext

Parameters indicated by a superscribed ∗ and † are taken from Bursztyn et al. (2007) and Cochran and Gao (2013), respectively. The other parameters are tuned to obtain a good fit between the simulated IUP and the clinical data in Buhimschi et al. (1997)

response. In this section, we investigate the impact of some of these parameters.

4.2.1 Excitation: spiking versus plateau APs The electrophysiological properties of USMCs undergo changes during pregnancy. These changes consist primarily of a progressive alteration in spike configuration from a single, plateau-type spike at early and mid-pregnancy, to a burst of APs (repetitive train) at term (Anderson et al. 1981; Kawarabayashi and Marshall 1981). The frequency of APs within a burst ranges from 0.1 to 10 Hz (Rabotti et al. 2010). In this simulation, two bursts of APs with higher frequency (0.35 Hz) at the beginning and lower (0.15 Hz) toward the end (Parkington et al. 1999), and with different durations are generated from the pacemakers’ region (upper-right side). The resulting IUPs are shown in Fig. 9. It is observed that by using the baseline diffusion coefficient in Table 2, u (and Vm ) will not propagate through the uterus. To overcome the problem, is increased from 0.6 × 10−3 to 0.1 m2 /s, and as a consequence, the propagation speed becomes ≈18 cm/s which is outside of the clinical range (Rabotti et al. 2010; Lammers 2013).

123

Fig. 9 The IUP generated based on spiking APs and the baseline IUP (solid line). The baseline is obtained from a plateau-type AP with depolarization time 90 s, while the spiking APs duration is 90 s (top insert) and 130 s (bottom insert)

4.2.2 Activation: Ca2+ desensitization Simultaneous measurements during spontaneous and electrically induced contractions of nonpregnant human myometria has shown that the force develops at a slower rate than [Ca2+ ]i increases and myosin light chains (MLC) are phosphorylated

Simulating uterine contraction by using an electro-chemo-mechanical model

Fig. 10 The effects of modifications in the myosin phosphorylation model (Sect. 2.2) on IUP. In the baseline k2 = 1.2387 s−1 has been utilized

(Word et al. 1994). Maximal and steady-state values of MLC phosphorylation were reached prior to the values of [Ca2+ ]i (Word et al. 1994) suggested that a desensitization of MLCK causes a decrease in the rate of MLC phosphorylation and force production. This can be simulated in two ways: by lowering the phosphorylation rate through a reduction in k1 (and k6 ), or by raising the dephosphorylation rate through an increase in k2 (and k5 ). We choose to increase k2 to twice its baseline value since it is simplest to implement. The result is shown in Fig. 10 together with the results for k2 equal to the baseline value and half the baseline value. The peak IUP decreases with increasing k2 , and vice versa, while the duration remains unchanged.

pregnant uteri (Sobotta 1891; Kawarabayashi and Marshall 1981; Weiss et al. 2006). To the best of authors’ knowledge, there is no information regarding to the fasciculi alignment in a pregnant uterus. This is unfortunate since any modification in these alignments will affect both the mechanical properties of the uterus and the electrophysiological activities of myometrium, including AP propagation and the [Ca2+ ]i dynamics. To show this, we simulate two different cases. First, the two USM families are rotated 5◦ while maintaining the right angle between them. Second, the first family is rotated 5◦ , while the second family is rotated −5◦ . The IUP for the first case does not differ from the baseline IUP, while for the second case, the IUP duration and peak value increase slightly (Fig. 11a). This may be caused by the longer depolarization duration (Fig. 11b) and/or the structural modification. 4.2.4 IUP sensitivity to κvol In addition to the physiologically related parameters, the effect of κvol in Eq. (18) on the IUP is also investigated herein. This investigation is motivated by noticing that the additive split of the Helmholtz free-energy function for an anisotropic material can lead to an unphysical behavior and affect the stress distribution considerably (Helfenstein et al. 2010; Annaidh et al. 2013; Weickenmeier and Jabareen 2014). Figure 12 illustrates how the IUP changes in response to different values of κvol .

5 Discussion 4.2.3 Contraction: the USM fasciculi alignment The USM fasciculi are assumed to be aligned longitudinally and circumferentially (Fig. 4) based on data from the structures of USM fasciculi and collagen fibers in non-

The need for a predictive tool to simulate uterine contractions motivated the work presented herein. An ECM model was developed to replicate the excitation, activation, and contraction of USMCs. The model is implemented in a 3D FE setting

Fig. 11 a The IUPs for different arrangements of USM fasciculi. b The AP for different arrangements of USM fasciculi evaluated at a point at the bottom of the domain

123

B. Sharifimajd et al.

Fig. 12 The IUP obtained for different values of κvol

to simulate uterine contractions during labor. The intrauterine pressure generated during the uterine contraction together with the uterine electrical activity were determined through the simulation and compared with experimental data. A parameter study was also carried out to investigate the impact of some important parameters on the model response. A homogenized structure is assumed for the uterine wall where the extracellular matrix deforms together with the smooth muscle. This implies that the extracellular material and the muscle cells present at a point X in the reference configuration must be present at the point x = x(X, t) in the current configuration, i.e., there is no relative motion between the constituents. Uterine smooth muscle fasciculi are interlaced with the extracellular matrix and motion between the constituents may occur. This coupling would contribute to the shear stress in the uterine wall, but this effect is not accounted for in the present study (see Holzapfel and Ogden 2009 for ideas in this direction). The reference configuration is assumed to be a truncated ellipsoid with dimensions approximately equal to that of the human uterus at 36 weeks of gestation. In this state, the material is assumed to be free of deformation. This approach neglects the residual strain introduced when deforming the uterus from its true stress-free state to the reference configuration. The residual strain is an important factor for the smooth muscle response since the generated muscle stress is length dependent. To compensate for the lack of residual opt strain in the USMCs, the parameter ε¯ a in Eq. (28) is set to a value >0 such that the maximum active force is generated for stretches well above those in the reference configuraopt tion. We have been unable to find a published value for ε¯ a in the uterus. It is, instead, computed using data for swine carotid artery reported in Murtada et al. (2012) and assuming ε¯ e = 0.05 at maximum smooth muscle force. The residual strain also has a significant impact on the stress field in the passive tissue and the question how to describe it arises. The residual strain in arteries is usually esti-

123

mated by measuring the opening angle obtained after cutting an arterial ring (Holzapfel et al. 2000). This stress-relieving cut is not easily transferred to the uterus which has a much more complex geometry and residual stress distribution. A good parametrization of the stress-free uterine state is yet to be published to the best of the authors’ knowledge. As a consequence, we choose to neglect residual strains. It is assumed that the strain-energy function Ψ is additively decoupled into a volumetric and a deviatoric part. This decomposition follows the tradition in nonlinear FE for nearly incompressible materials. Serious critique regarding its validity for anisotropic materials has been raised (Annaidh et al. 2013). An example of the problems that may be encountered when using an additive decomposition is given in Helfenstein et al. (2010). They showed that growing lateral stretches cause a stress reduction when extending a cube along the reinforcing fibers. The unrealistic response is manifested by a volume ratio J > 1 and is intimately linked to the parameter κvol . It is therefore important to adjust this parameters such that the material behaves (nearly) incompressible. The baseline simulation described in Sect. 4.1 was carried out using κvol = 1 MPa which gave a peak J of 1.03. It was later observed that a value around 10 MPa, or higher, might have been more appropriate, assuming the uterine wall to be perfectly incompressible. For this value J was very close to one. The outer boundary of the uterus is subjected to a uniform pressure taken to be 28 mmHg. This is a gross simplification of the real situation where the surrounding abdominal tissue will cause a complex and time-varying load distribution. A uniform pressure is the simplest possible nonzero boundary condition, however, and it introduces the physiological uterine tone. Furthermore, the cervical end of the uterus is assumed to be fixed in all directions. In reality the cervix is a compliant structure (Badir et al. 2013; Mazza et al. 2014; Briggs et al. 2015; Myers et al. 2015), and a more realistic, but still simple, boundary condition would be to assume a linearly elastic support. Numerical experiments using data from Badir et al. (2013, Fig 11) indicate that the IUP remains relatively unchanged by this condition, but unphysiological oscillations in the displacement of the uterine wall led to the conclusion that this approach needs further work. Future work will include further validation of the model against measurements of important clinical factors such as IUP. More experimental data would be useful to restrict the possible ranges of parameters in the model. This include geometrical data such as USM and collagen fibers architecture, which can be obtained noninvasively (Weiss et al. 2006), and mechanical properties of the cervix (Badir et al. 2013; Mazza et al. 2014) for implementation of more realistic boundary conditions. A detailed statistical sensitivity analysis for the simulation is also desirable, possibly using response surface

Simulating uterine contraction by using an electro-chemo-mechanical model

methodology (Myers et al. 2009) as a way of reducing the amount of time spent on FE simulations. Acknowledgments The authors would like to thank the anonymous reviewers for their valuable comments which helped to improve the paper. We also thank Daniel Ericsson and other members of the support team at Comsol Multiphysics for their guidance and Prof. Anders Klarbring for the helpful discussion.

References Aliev RR, Panfilov AV (1996) A simple two-variable model of cardiac excitation. Chaos Solitons Fractals 7(3):293–301 Ambrosi D, Arioli G, Nobile F, Quarteroni A (2011) Electromechanical coupling in cardiac dynamics: the active strain approach. SIAM J Appl Math 71(2):605–621 Anderson GF, Kawarabayashi T, Marshall JM (1981) Effect of indomethacin and aspirin on uterine activity in pregnant rats: comparison of circular and longitudinal muscle. Biol Reprod 24(2):359–372 Annaidh AN, Destrade M, Gilchrist M, Murphy J (2013) Deficiencies in numerical models of anisotropic nonlinearly elastic materials. Biomech Model Mechanobiol 12(4):781–791 Arner A (1982) Mechanical characteristics of chemically skinned guinea-pig taenia coli. Pflügers Arch 395(4):277–284 Badir S, Bajka M, Mazza E (2013) A novel procedure for the mechanical characterization of the uterine cervix during pregnancy. J Mech Behav Biomed Mater 27:143–153 Bastos LF, Lobo MF, van Meurs WL, Ayres-de Campos D (2010) An intrauterine pressure generator for educational simulation of labour and delivery. Med Eng Phys 32(7):740–745 Benson AP, Clayton RH, Holden AV, Kharche S, Tong WC (2006) Endogenous driving and synchronization in cardiac and uterine virtual tissues: bifurcations and local coupling. Philos Trans R Soc A Math Phys Eng Sci 364(1842):1313–1327 Blencowe H, Cousens S, Oestergaard MZ, Chou D, Moller AB, Narwal R, Alder A, Vera Garcia C, Rohde S, Say L, Lawn JE (2012) National, regional, and worldwide estimates of preterm birth rates in the year 2010 with time trends since 1990 for selected countries: a systematic analysis and implications. Lancet 379(9832):2162– 2172 Böl M, Schmitz A, Nowak G, Siebert T (2012) A three-dimensional chemo-mechanical continuum model for smooth muscle contraction. J Mech Behav Biomed Mater 13:215–229 Briggs BN, Stender ME, Muljadi PM, Donnelly MA, Winn VD, Ferguson VL (2015) A hertzian contact mechanics based formulation to improve ultrasound elastography assessment of uterine cervical tissue stiffness. J Biomech 48(9):1524–1532 Bugg GJ, Stanley E, Baker PN, Taggart MJ, Johnston TA (2006) Outcomes of labours augmented with oxytocin. Eur J Obstet Gynecol Reprod Biol 124(1):37–41 Buhimschi C, Boyle MB, Garfield RE (1997) Electrical activity of the human uterus during pregnancy as recorded from the abdominal surface. Obstet Gynecol 90(1):102–111 Buhimschi CS, Buhimschi IA, Malinow AM, Weiner CP (2003) Myometrial thickness during human labor and immediately post partum. Am J Obstet Gynecol 188(2):553–559 Buhimschi CS, Buhimschi IA, Malinow AM, Weiner CP (2004) Intrauterine pressure during the second stage of labor in obese women. Obstet Gynecol 103(2):225–230 Burdyga T, Wray S, Noble K (2007) In situ calcium signaling. Ann NY Acad Sci 1101(1):85–96

Bursztyn L, Eytan O, Jaffa AJ, Elad D (2007) Mathematical model of excitation-contraction in a uterine smooth muscle cell. Am J Physiol Cell Physiol 292(5):C1816–C1829 Cochran AL, Gao Y (2013) A model and simulation of uterine contractions. Math Mech Solids 20:540–564 Csapo A, Pinto-Dantas C (1965) The cyclic activity of the nonpregnant human uterus. A new method for recording intrauterine pressure. Fertil Steril 17(1):34–38 Danforth DN (1971) Textbook of obstetrics and gynecology, by 54 authors. Medical Book Dept., Harper & Row FitzHugh R (1961) Impulses and physiological states in theoretical models of nerve membrane. Biophys J 1(6):445–466 Gifford DS, Morton SC, Fiske M, Keesey J, Keeler E, Kahn KL (2000) Lack of progress in labor as a reason for cesarean. Obstet Gynecol 95(4):589–595 Göktepe S, Kuhl E (2010) Electromechanics of the heart: a unified approach to the strongly coupled excitation–contraction problem. Comput Mech 45(2–3):227–243 Hai CM, Murphy RA (1988) Cross-bridge phosphorylation and regulation of latch state in smooth muscle. Am J Physiol Cell Physiol 254(1):C99–C106 Helfenstein J, Jabareen M, Mazza E, Govindjee S (2010) On nonphysical response in models for fiber-reinforced hyperelastic materials. Int J Solids Struct 47(16):2056–2061 Herlihy JT, Murphy RA (1973) Length–tension relationship of smooth muscle of the hog carotid artery. Circ Res 33(3):275–283 Hodgkin AL, Huxley AF (1952) A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 117(4):500 Holzapfel GA (2000) Nonlinear solid mechanics: a continuum approach for engineering. Wiley, West Sussex Holzapfel GA, Ogden RW (2009) Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philos Trans R Soc Lond A Math Phys Eng Sci 367(1902):3445–3475 Holzapfel GA, Gasser TC, Ogden RW (2000) A new constitutive framework for arterial wall mechanics and a comparative study of material models. J Elast Phys Sci Solids 61(1–3):1–48 Huxley A (1974) Muscular contraction. J Physiol 243(1):1–43 Johnson T, Clayton C (1955) Studies in placental action during prolonged and dysfunctional labours using radioactive sodium. BJOG Int J Obstet Gynaecol 62(4):513–522 Kawarabayashi T, Marshall JM (1981) Factors influencing circular muscle activity in the pregnant rat uterus. Biol Reprod 24(2):373–379 Keener JP, Sneyd J (2008) Mathematical physiology: I: cellular physiology, vol 1. Springer, New York Korogod S, Kochenov A (2013) Mathematical model of the calciumdependent chloride current in a smooth muscle cell. Neurophysiology 45(5–6):369–378 La Rosa PS, Eswaran H, Preissl H, Nehorai A (2012) Multiscale forward electromagnetic model of uterine contractions during pregnancy. BMC Med Phys 12(1):4 Lammers WJ (2013) The electrical activities of the uterus during pregnancy. Reprod Sci 20(2):182–189 Londino S (2014) Birth root midwifery. http://www.sandralondino. com/stages-of-labor/ Maul H, Maner WL, Saade GR, Garfield RE (2003) The physiology of uterine contractions. Clin Perinatol 30(4):665–676 Mazza E, Parra-Saavedra M, Bajka M, Gratacos E, Nicolaides K, Deprest J (2014) In vivo assessment of the biomechanical properties of the uterine cervix in pregnancy. Prenat Diagn 34(1):33–41 McMahon TA (1984) Muscles, reflexes, and locomotion. Princeton University Press, Princeton Murtada SI, Holzapfel GA (2014) Investigating the role of smooth muscle cells in large elastic arteries: a finite element analysis. J Theor Biol 358:1–10

123

B. Sharifimajd et al. Murtada SI, Kroon M, Holzapfel GA (2010) A calcium-driven mechanochemical model for prediction of force generation in smooth muscle. Biomech Model Mechanobiol 9(6):749–762 Murtada SI, Arner A, Holzapfel GA (2012) Experiments and mechanochemical modeling of smooth muscle contraction: significance of filament overlap. J Theor Biol 297:176–186 Mwaniki MK, Atieno M, Lawn JE, Newton CR (2012) Long-term neurodevelopmental outcomes after intrauterine and neonatal insults: a systematic review. Lancet 379(9814):445–452 Myers KM, Feltovich H, Mazza E, Vink J, Bajka M, Wapner RJ, Hall TJ, House M (2015) The mechanical role of the cervix in pregnancy. J Biomech 48(9):1511–1523 Myers RH, Montgomery DC, Anderson-Cook CM (2009) Response surface methodology: process and product optimization using designed experiments, vol 705. Wiley, New York Nagumo J, Arimoto S, Yoshizawa S (1962) An active pulse transmission line simulating nerve axon. Proc IRE 50(10):2061–2070 Nobile F, Quarteroni A, Ruiz-Baier R (2012) An active strain electromechanical model for cardiac tissue. Int J Numer Methods Biomed Eng 28(1):52–71 Padubidri V, Anand E (2006) Textbook of obstetrics. BI Publications Pvt Ltd., New Delhi Parkington HC, Tonta MA, Davies NK, Brennecke S, Coleman H (1999) Hyperpolarization and slowing of the rate of contraction in human uterus in pregnancy by prostaglandins e2 and f2α: involvement of the na+ pump. J Physiol 514(1):229–243 Parthimos D, Edwards D, Griffith T (1999) Minimal model of arterial chaos generated by coupled intracellular and membrane ca2+ oscillators. Am J Physiol Heart Circ Physiol 277(3):H1119–H1144 Rabotti C, Mischi M, Oei SG, Bergmans JW (2010) Noninvasive estimation of the electrohysterographic action-potential conduction velocity. IEEE Trans Biomed Eng 57(9):2178–2187 Rumpel T, Schweizerhof K, Haßler M (2005) Efficient finite element modelling and simulation of gas and fluid supported membrane and shell structures. Textile composites and inflatable structures. Springer, Netherlands, pp 153–172 Sandow A (1952) Excitation-contraction coupling in muscular response. Yale J Biol Med 25(3):176 Sharifimajd B, Stålhand J (2013) A continuum model for skeletal muscle contraction at homogeneous finite deformations. Biomech Model Mechanobiol 12(5):965–973 Sharifimajd B, Stålhand J (2014) A continuum model for excitationcontraction of smooth muscle under finite deformations. J Theor Biol 355:1–9 Simo JC, Taylor RL (1991) Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms. Comput Methods Appl Mech Eng 85(3):273–310 Snell RS (2011) Clinical anatomy by regions. Lippincott Williams & Wilkins, Wolters Kluwer Sobotta J (1891) Beiträge zur vergleichenden anatomie und entwickelungsgeschichte der uterusmuskulatur. Archiv für Mikroskopische Anatomie 38(1):52–100

123

Spencer AJM et al (1984) Continuum theory of the mechanics of fibrereinforced composites, vol 282. Springer, New York Stålhand J, Klarbring A, Holzapfel GA (2011) A mechanochemical 3d continuum model for smooth muscle contraction under finite strains. J Theor Biol 268(1):120–130 Tong WC, Choi CY, Karche S, Holden AV, Zhang H, Taggart MJ (2011) A computational model of the ionic currents, ca2+ dynamics and action potentials underlying contraction of isolated uterine smooth muscle. PloS One 6(4):e18,685 Toth TI, Schmidt J, Büschges A, Daun-Gruhn S (2013) A neuromechanical model of a single leg joint highlighting the basic physiological role of fast and slow muscle fibres of an insect muscle system. PloS One 8(11):e78,247 Usher R, McLean F (1969) Intrauterine growth of live-born caucasian infants at sea level: standards obtained from measurements in 7 dimensions of infants born between 25 and 44 weeks. J Pediatr 74(6):901–910 Uvelius B (1976) Isometric and isotonic length-tension relations and variations in cell length in longitudinal smooth muscle from rabbit urinary bladder. Acta Physiol Scand 97(1):1–12 Victorri B, Vinet A, Roberge FA, Drouhard JP (1985) Numerical integration in the reconstruction of cardiac action potentials using hodgkin-huxley-type models. Comput Biomed Res 18(1):10–23 Weickenmeier J, Jabareen M (2014) Elastic-viscoplastic modeling of soft biological tissues using a mixed finite element formulation based on the relative deformation gradient. Int J Numer Methods Biomed Eng 30(11):1238–1262 Weiss S, Jaermann T, Schmid P, Staempfli P, Boesiger P, Niederer P, Caduff R, 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 Cell Evolut Biol 288(1):84–90 Word RA, Tang DC, Kamm KE (1994) Activation properties of myosin light chain kinase during contraction/relaxation cycles of tonic and phasic smooth muscles. J Biol Chem 269(34):21,596–21,602 Yang J, Clark JW Jr, Bryan RM, Robertson C (2003) The myogenic response in isolated rat cerebrovascular arteries: smooth muscle cell model. Med Eng Phys 25(8):691–709 Young RC, Hession RO (1999) Three-dimensional structure of the smooth muscle in the term-pregnant human uterus. Obstet Gynecol 93(1):94–99 Zhu BP, Grigorescu V, Le T, Lin M, Copeland G, Barone M, Turabelidze G (2006) Labor dystocia and its association with interpregnancy interval. Am J Obstet Gynecol 195(1):121–128